Linear analysis on the onset of thermal convection of highly compressible fluids with variable viscosity and thermal conductivity in spherical geometry: implications for the mantle convection of super-Earths

In this paper, we carried out a series of linear analyses on the onset of thermal convection of highly compressible fluids whose physical properties strongly vary in space in convecting vessels either of a three-dimensional spherical shell or a two-dimensional spherical annulus geometry. The variations in thermodynamic properties (thermal expansivity and reference density) with depth are taken to be relevant for the super-Earths with ten times the Earth’s mass, while the thermal conductivity and viscosity are assumed to exponentially depend on depth and temperature, respectively. Our analysis showed that, for the cases with strong temperature dependence in viscosity and strong depth dependence in thermal conductivity, the critical Rayleigh number is on the order of 108–109, implying that the mantle convection of massive super-Earths is most likely to fall in the stagnant-lid regime very close to the critical condition, if the properties of their mantle materials are quite similar to the Earth’s. Our analysis also demonstrated that the structures of incipient flows of stagnant-lid convection in the presence of strong adiabatic compression are significantly affected by the depth dependence in thermal conductivity and the geometries of convecting vessels, through the changes in the static stability of thermal stratification of the reference state. When the increase in thermal conductivity with depth is sufficiently large, the thermal stratification can be greatly stabilized at depth, further inducing regions of insignificant fluid motions above the bottom hot boundaries in addition to the stagnant lids along the top cold surfaces. We can therefore speculate that the stagnant-lid convection in the mantles of massive super-Earths is accompanied by another motionless regions at the base of the mantles if the thermal conductivity strongly increases with depth (or pressure), even though their occurrence is hindered by the effects the spherical geometries of convecting vessels.


Introduction
Super-Earths are extra-solar planets which have small masses (compared to that of Jupiter) and high mean density larger than 5000 kg/m 3 (Howard et al. 2010). Since the first discovery of such massive terrestrial planets (Rivera et al. 2005), a growing effort has been devoted to the studies in the internal (radial) structure (e.g., Seager Kameyama Earth, Planets and Space (2021) 73:167 et al. 2007;Valencia et al. 2007b;Sotin et al. 2007;Wagner et al. 2011), the thermal evolution (e.g., Papuc and Davies 2008;Kite et al. 2009;Korenaga 2010;Gaidos et al. 2010;Tachinami et al. 2011;Čížková et al. 2017) and the habitability (e.g., Noack et al. 2017;Tosi et al. 2017;Foley 2019) of super-Earths. During the last decades, the studies on mantle dynamics have been playing an important role in addressing these issues (e.g., Shahnas et al. 2018;van den Berg et al. 2018).
It can be intuitively understood that, because of their large sizes, the convection in the mantles of massive super-Earths is most likely to occur in a quite different manner from those in the Earth and smaller terrestrial planets. For example, the variations of physical properties of mantle materials are expected to be larger in the much thicker mantles of super-Earths than in the Earth's (e.g., Karato 2011). Indeed, the effects of strong spatial variations in physical properties (such as viscosity and thermal conductivity) on the convecting flow patterns have been widely studied so far, under the assumption that the effects of other physical mechanisms such as compressibility can be minor (e.g., Kameyama and Ogawa 2000;Miyauchi and Kameyama 2013). In particular, the effects of strong temperature dependence in viscosity have been intensively studied by numerical models of convection in three-dimensional spherical geometry (e.g., Yao et al. 2014;Yanagisawa et al. 2016;Guerrero et al. 2018).
On the other hand, it has been well acknowledged that the effect of adiabatic compression affects the nature of thermal convection (Jarvis and McKenzie 1980;King et al. 2010;Liu and Zhong 2013). This effect is expected to be of great importance in the convection in the thick mantles of massive super-Earths, considering that the pressure can reach on the order of terapascals in their interiors (Tsuchiya and Tsuchiya 2011). Indeed, our earlier studies (Kameyama and Kinoshita 2013;Kameyama et al. 2015;Kameyama 2016) showed that the nature of thermal convection in the thick mantles of massive super-Earths is strongly affected by the interplay between the effects of adiabatic compression and those of spatial variations in physical properties of mantle materials. However, in these studies, the effects of spherical geometries of planetary mantles have been thoroughly ignored, since the analyses had been carried out using the convecting vessels with the Cartesian (planar) domains. It is therefore very important to estimate the effects of the strong adiabatic compression and those of strong variations in physical properties of mantle materials on the convecting flow structures in spherical geometries, in order to gain the fundamental insights into the mantle convection of massive super-Earths.
In this study, we carry out a series of linear analyses on the thermal convection in spherical domains of a highly compressible fluid with the spatial variation in physical properties, to investigate how the spherical geometries of planetary mantles can affect the critical state and condition of thermal convection in the presence of the interplay between the adiabatic compression and spatial variations in physical properties such as viscosity and thermal conductivity. In particular, we will focus on the influences of adiabatic changes in temperature on the incipient convection, by using a simple 1D model of a basally heated fluid. By further comparing the present results with our earlier one in the Cartesian (planar) model (Kameyama 2016), we aim at deepening the insights into the thermal convection in the mantles of massive terrestrial super-Earths both from the thermodynamic and fluid-dynamic aspects.

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

Physical properties of modeled compressible fluid
As a model of pressure dependence in thermodynamic properties of mantle materials of massive super-Earths, we employ the same variations in thermal expansivity α and reference density ρ with depths as in our earlier works Miyagoshi et al. 2014Miyagoshi et al. , 2015Kameyama 2016;Kameyama and Yamamoto 2018). The thermal expansivity α decreases with depth, and its radial profile is taken to be similar to that estimated for MgO from an ab initio calculation (Tsuchiya and Tsuchiya 2011) for the mantles of super-Earths with ten times the Earth's mass. That is, the value of α(r = r min ) ≡ α bot at the bottom boundary is as small as about 1/14 of the value α(r = r max ) ≡ α top ( = 4 × 10 −5 K −1 ) at the top surface, and the value of its average α ave is taken to be On the other hand, the reference density ρ is assumed to linearly increase with depth, yielding the value of ρ(r = r min ) at the bottom boundary larger than ρ(r = r max ) ≡ ρ top at the top surface by about 3.17 (Valencia et al. 2006). We note that the linear increase in ρ assumed here is not likely to spoil the significance of the present analysis despite its geophysical unsoundness, since the vertical profile of (reference) density has a very minor effect on the nature of critical states of thermal convection of highly compressible fluids .
In addition to the variations in thermodynamic properties ( α and ρ ), we also assume the spatial variations in transport properties. In this study, the thermal conductivity k is taken to be exponentially dependent on depth as, where k top is the thermal conductivity at r = r max , and r k is a dimensionless constant equal to the ratio of k(r = r min ) ≡ k bot to k top . For a geophysical significance, we restrict ourselves to the cases where r k ≥ 1 ; the thermal conductivity k increases with depth for r k > 1 . According to the recent estimates both by theoretical calculations (Dekura et al. 2013;Tsuchiya 2017, 2019) and by laboratory experiments (Haigis et al. 2012;Ohta et al. 2012Ohta et al. , 2017Lobanov et al. 2017), the thermal conductivity of mantle materials can increase by up to about one-two orders of magnitude across the mantles of super-Earths. In this study, we varied r k up to 100.
On the other hand, the viscosity η is assumed to be exponentially dependent on temperature T as, where η bot is the viscosity for T = T bot , and r η is a dimensionless constant given by r η = η(T = T top )/η(T = T bot ) . In (3), for simplicity, we ignored the dependence on pressure (or depth) in viscosity, since it is still highly controversial if the viscosity of mantle materials increases or not as the pressure increases up to the levels relevant to the interior of super-Earths (Stamenković et al. 2011;Karato 2011;Wagner et al. 2012;Tackley et al. 2013;Noack and Breuer 2014). In this study, we restrict ourselves to the cases where r η ≥ 1 ; the viscosity η decreases with temperature for r η > 1 . Using this assumption, we will focus on the effects of the reduction of viscosity in the deep mantle of the super-Earths proposed in several earlier studies (Karato 2011;Stein et al. 2011;Wagner et al. 2012;Tackley et al. 2013).
Here, some comments are necessary on the depth dependence in thermal conductivity k in this study. First, the exponential increase in k with depth given by (2) is chosen simply from analytical and numerical convenience as in our earlier studies (Kameyama and Kinoshita 2013;Kameyama 2016;Kameyama and Yamamoto 2018) and, in other words, may not be necessarily supported by studies on mineral physics conducted so far (e.g., Dalton et al. 2013;Hsieh et al. 2017Hsieh et al. , 2018Deschamps and Hsieh 2019). Second, it is well acknowledged that the thermal conductivity k depends not only on depth (or pressure) but also on other parameters such as composition and temperature. For example, the thermal conductivity k of ferro-periclase is much higher and more sensitive to pressure (about a factor of 10 increase in the pressure range of the Earth's mantle) than that of bridgmanite. Similarly, the values of k of both minerals are sensitive to the iron content. In addition, the temperature dependence in k is expected to compensate for the increase in k with depth, since high temperatures tend to reduce the contributions of heat transport by lattices (or phonons). In short, the present study may overestimate the values of r k ≡ k bot /k top in silicate mantles of super-Earths by ignoring the dependence of k on composition and temperature.

Basic equations
Under the truncated anelastic liquid approximation (TALA), the basic equations can be written in a dimensionless form as: where the quantities with overbars denote those of reference state which vary only in the radial direction, and the nondimensionalization is carried out with the length scale of d ≡ r max − r min , time scale of d 2 /κ top (where κ top = k top /ρ top C p is the value of thermal diffusivity at the top surface), and temperature scale of T ≡ T bot − T top . Here, ⊗ stands for the tensor product of vectors, v ⊗ ∇ is the transpose of ∇ ⊗ v , I is the identity tensor, e r is the unit vector in radial direction positive outward, v is the velocity vector, v r ≡ (v · e r ) is the velocity in the radial direction, p is pressure, and ε II is the second invariant of strain rate tensor. In the right-hand side of (6), the first, second, and third terms represent the (4) 0 = ∇ · (ρv), effects of apparent heat transport (advection plus conduction), adiabatic temperature change, and viscous dissipation (frictional heating), respectively. The present formulation includes two nondimensional parameters. The first one is the dissipation number Di of the modeled fluid defined by, The second one is the Rayleigh number Ra of thermal convection defined by, Note here that, as indicated by (3), the viscosity is scaled by the value η bot at the bottom boundary while for other properties the reference values are those at the top surface. The values of Di and Ra for the mantles of terrestrial planets are expected to increase with the mass M of the planets, coming from the increase in the mass, gravity, and the mantle thickness of the terrestrial planets whose composition is similar to the Earth's (e.g., Seager et al. 2007;Valencia et al. 2007a;Wagner et al. 2011). According to the scaling relationship [ g ∝ M 0.5 and d ∝ M 0.28 ; Valencia et al. (2007a)] together with Di ∼ 0.7 and Ra = O(10 7 ) for the Earth's mantle (Schubert et al. 2001), one can estimate Di ∼ 5 and Ra = O (10 8 -10 9 ) for the mantles of super-Earths whose mass is ten times larger than the Earth's Miyagoshi et al. 2014Miyagoshi et al. , 2015Miyagoshi et al. , 2018. In the present study, we fixed Di = 5 , while the values of Ra estimated above can be compared with those of the critical Rayleigh numbers derived from the linear stability analysis.
In the linear stability analysis presented in this paper, we solve for the temporal evolution of an infinitesimal perturbation superimposed to a reference state described by a stationary (motionless) state with steady one-dimensional heat conduction in the radial direction. In the following, the quantities with overbars and primes denote those of reference state and perturbation, respectively. The dimensionless equations for the reference state are written as, while those for infinitesimal perturbations are given by, In deriving these equations, the second-order terms of infinitesimal (primed) quantities are ignored. This is the reason why the effect of viscous dissipation is eliminated in (14).
At the top and bottom boundaries, we kept the perturbations in temperature T ′ and the radial velocity v ′ r to be zero. For the horizontal velocity at the boundaries, on the other hand, we assumed a free-slip condition.

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

Static stability of thermal stratification for reference (steady) states
We first see how the structure and static stability of the thermal stratification of the reference state vary depending on the model geometries and on the depth dependence in thermal conductivity k . By the solid lines in Fig. 1a, we show the radial profiles of reference temperature T (r) derived from Eq. (10) for different model geometries and various values of r k . [Note that the temperature dependence in viscosity η affects neither the thermal state nor its static stability of the reference state given by Eq. (10)]. Figure 1b shows the magnitude of its radial gradient |dT /dr| . To further highlight the effects of model geometries, we also show in Fig. 1 the similar plots of T and its vertical gradient dT /dz for the Cartesian geometry derived in our earlier work (Kameyama 2016). From the plots in Fig. 1a, b, we can see that the radial profiles of reference temperature T are significantly affected by the values of r k as well as by the model geometries. The comparison of the plots for different r k and given geometries shows that the values of T in the fluid layer are higher for larger r k . The increase in T for large r k comes from the changes in the radial profiles of |dT /dr| within the fluid layer depending on r k [see also Eq. (10)]: a larger increase in k with depth induces a larger decrease in |dT /dr| with depth. This in turn makes |dT /dr| larger in the shallow cold part and smaller in the deep hot part of the fluid layer (see Fig. 1b), which consequently raises T in the entire layer (see Fig. 1a). The comparison of the plots of T for 3D and 2D spherical geometries, on the other hand, shows that for a given r k the values of T in the fluid layer are lower in 3D geometry than in 2D geometry. The difference in T due to the model geometries comes from the difference in the changes in the surface area with depth. The decrease in the surface area with depth is more significant in 3D cases ( ∝ r −2 ) than in 2D cases ( ∝ r −1 ). Owing to this effect, the magnitude of |dT /dr| in the shallow part of fluid layer is smaller in 3D cases than in 2D cases and, at the same time, the magnitude in the deep part is larger (see Fig. 1b). This results in lower T in the entire layer for 3D cases than in 2D cases. The effect of the changes in the surface area with depth can be further confirmed from the comparison with the Cartesian cases (Kameyama 2016). Among the three model geometries presented in Fig. 1, indeed, the values of T for given r k are the highest for Cartesian cases where the effect of the changes in the surface area with depth is absent.
To see the static stability of the thermal stratification more clearly, we show in Fig. 1c the radial profiles of potential temperature T pot of the reference states. Here, T pot at the dimensionless position r is calculated by We also show by thick gray lines the ranges of depths where T pot increases with height. Note that in these regions the thermal stratification given by T is stable, as is inferred from the Brunt-Väisälä frequency N: Here the dimensionless value of N is given in the present numerical configuration by where a nondimensionalization is carried out by a frequency scale of g/d , and (dT /dr) s is the adiabatic temperature gradient whose dimensionless value is defined by In the regions where N 2 > 0 (i.e., N is real), the thermal stratification is stable against an imposed infinitesimal displacement of fluid parcels in the direction of gravity. Indeed, by differentiating (15) with respect to r we get which indicates that dT pot /dr is positive when |dT /dr| s > |dT /dr| and, hence, N 2 > 0 . In addition, the static stability of thermal stratification is stronger for larger N 2 and dT pot /dr.
The comparison between the plots in Fig. 1c clearly shows that the static stability of the thermal stratification is affected not only by the values of r k but also by the geometries of fluid layer. For both 3D and 2D spherical geometries, regions of stable thermal stratification with dT pot /dr > 0 emerge in the deeper part of the fluid layer when r k ≥ 30 . This is because a large r k reduces the magnitude of |dT /dr| in the deeper part and, hence, enhances the static stability of thermal stratification [see also Eq. (18)]. These regions become thicker and the thermal stratification there becomes more stable for larger r k . We note, on the other hand, from the comparison between 3D and 2D cases that for given r k ≥ 30 the regions of stable stratification are thinner for 3D cases than those for 2D cases. This is again due to the difference in the changes in the surface area with depth. For 3D cases, the magnitude of |dT /dr| in the deeper part of the fluid layer becomes larger than in 2D cases, which in turn reduces the static stability of thermal stratification at depths. The effect of model geometries on the static stability can be further confirmed from the comparison with the Cartesian cases (Kameyama 2016) where the emergence of regions of stable stratification is greatly enhanced with respect to that with spherical cases.
We will conclude this subsection by estimating the condition under which the regions of stable stratification Plots against height of the magnitude of the gradient of temperature |dT /dr| (or |dT /dz| ) in the radial (vertical) directions. c Plots against height of the potential temperature T pot . Also shown by thick gray lines are the ranges of depths where the thermal stratification given by T is stable owing to the increase in T pot with height emerge at the bottom of the fluid layer. This condition is given by dT /dr s > dT /dr at r = r min . Using (17), this inequality can be further rewritten as In Table 1, we summarize the values of dT /dr r=r min calculated for several r k and 3D spherical shell geometry with r min /r max = 0.5 (see Fig. 1a, b), together with the resulting values of Di thres . In calculating Di thres we used α bot ∼ 1/14 and T bot = 1.1 in dimensionless units. In Table 1 we also show the values of the threshold mass M thres of planets above which Di > Di thres in their mantles, based on the scaling relationship of Di ∝ gd ∝ M 0.78 described in the previous section. Note, however, that in these estimates we employed the same decrease in thermal expansivity α with depth as that for M = 10M E (where M E is the Earth's mass). That is, these estimates may not be very accurate for the cases with M = 10M E . In particular, the values of Di thres and M thres are most likely to be smaller than those estimated here for the cases with M < 10M E , where both α ave /α top and α bot /α top are expected to be larger than those for M = 10M E and, in other words, the effects of adiabatic compression are more significant.

Linear stability analyses for the onset of thermal convection
Next we carry out linear stability analyses for the onset of thermal convection of highly compressible fluids in the presence of spatial variations in physical properties (temperature dependence of viscosity η and depth dependence in thermal conductivity k) in the convecting vessels with different geometries. In particular, we compared the changes in the critical states and structures of incipient flows with those in the static stability of thermal stratification for the reference state, to investigate how the nature of incipient convection of compressible fluids with spatial variations in physical properties is affected by the adiabatic temperature change and, in other words, the occurrence of regions of stable stratification. In Fig. 2, we show the plots against the viscosity contrast r η across the fluid layer of (a) the absolute minimum value of critical Rayleigh number Ra c0 and (b) the spherical harmonic degree ℓ c0 or the horizontal wavenumber K c0 corresponding to Ra c0 , obtained for different model geometries and various values of r k . Figure 2a shows that in general the values of Ra c0 become larger for larger r η and r k , which can be intuitively understood from an overall increase in viscosity and thermal conductivity with r η and/or r k . Figure 2b shows, on the other hand, that the values of ℓ c0 and K c0 rapidly increase with r η for r η O (10 3 -10 4 ), regardless of the model geometries and the values of r k . This increase in ℓ c0 or K c0 with r η has been recognized by the transition of "stagnant-lid" regime of thermal convection (Stengel et al. 1982;. Taken together with the plots in Fig. 2a, the values of Ra c0 for the stagnant-lid regime and r k ≥ 30 are as large as O (10 8 -10 9 ), which is very close to the estimate of the Rayleigh number of mantle convection of super-Earths whose mass is about ten times the Earth's. This suggests that in massive super-Earths the stagnant-lid type of mantle convection may occur under the conditions very close to the critical ones. In addition, among the three geometries presented here, the 3D spherical geometry needs the largest values of r η for the transition into "stagnant-lid" regime, which is consistent with an earlier finding by linear ) and nonlinear (Guerrero et al. 2018) studies on the effects of curvature of fluid layer.
We next focus on the nature of critical conditions for strong temperature dependence in viscosity ( r η O (10 3 -10 4 ) where the thermal convection falls in the stagnant-lid regime. For both the 2D and 3D spherical geometries, the changes in Ra c0 and those in ℓ c0 and K c0 with r η are different depending on whether r k ≤ 10 or r k ≥ 30 and, hence, whether the regions of stable thermal stratification occur or not in the reference state (see also Fig. 1). For the cases with r k ≤ 10 lacking in stable regions, the values of Ra c0 become almost independent of r k for large r η (Fig. 2a), and the values of ℓ c0 and K c0 are smaller for larger r k (Fig. 2b). For the cases with r k ≥ 30 in the presence of stable regions, by contrast, the values of Ra c0 are always larger for larger r k , and the values of ℓ c0 and K c0 are smaller for larger r k . We also note from the plots for the Cartesian geometry that the similar changes in Ra c0 and K c0 with r η occur depending on whether r k ≤ 10 or Table 1 The values of the temperature gradient |dT /dr| r=r min at the bottom of the fluid layer of 3D spherical shell of r min /r max = 0.5 calculated with several values of r k = k bot /k top , the threshold values of dissipation number Di thres and planetary mass M thres divided by the Earth's mass M E above which regions of stable thermal stratification occur at the bottom of the fluid layer  . 2 The plots against r η of a the absolute minimum value of critical Rayleigh number Ra c0 and b the spherical harmonic degree ℓ c0 or horizontal wavenumber K c0 of perturbation corresponding to Ra c0 , obtained for different model geometries and various values of r k in the range of 1 ≤ r k ≤ 100 r k ≥ 30 and, in other words, whether the regions of stable thermal stratification exist or not.
To see how the depth dependence in thermal conductivity k affects the critical states of stagnant-lid convection, we show in Fig. 3 the schematic illustration on two-dimensional planes of the incipient flows obtained for several values of r k and various model geometries together with r η = 10 8 . (For the cases in 3D spherical geometry, we show the flow structures in meridian planes assuming the incipient flows occurring with only zonal components whose spherical harmonic order is m = 0 .) Shown in color are the distributions in the perturbations in temperature T ′ whose color scales are indicated at the bottom of the figure. Shown by solid lines are the contours of a "potential" ψ of the mass flux ρv ′ which satisfies where e n is the unit vector perpendicular to the plane. (Note that the function ψ is a generalization of the stream function for the Boussinesq or extended Boussinesq approximation where ρ = 1 .) As can be seen from Fig. 3a, the incipient flows for r k = 1 occur as a wellknown stagnant-lid convection. Indeed, thick lids of cold and highly viscous fluid develop along the top surface where the fluid motion is insignificant. The convection occurs only beneath the thick and stagnant lids, which results in a small vertical extent of convection cells. The horizontal extent of convection cells also becomes much smaller than the thickness of entire layer, in accordance with the reduction in their vertical extent. This results in large ℓ c0 and K c0 of the infinitesimal perturbation (see also Fig. 2b).
The comparison of Fig. 3b and c with 3a clearly shows that the structures of incipient flows of stagnant-lid convection are significantly affected by the increase in r k , particularly through the changes in the flow strength not only in uppermost part but also in the lowermost part of fluid layer. On the one hand, the changes in the flow strength in the uppermost part is characterized by the thinning of stagnant lids along the top surface by the increase in r k . This is because for larger r k the temperature T of the reference state is higher (see also Fig. 1a) and, hence, the regions of high viscosity η become (20) ρv ′ = ∇ × (ψe n ), thinner owing to the strong temperature dependence in viscosity. On the other hand, the change in the lowermost part can be seen as an emergence of the regions of insignificant fluid motion just above the bottom boundary when r k is as large as 100. The emergence of such distinct regions for very large r k results from the changes in the static stability of thermal stratification for the reference state. As can be seen from Fig. 1c, thick regions of stable thermal stratification occur above the bottom boundaries for r k = 100 , while no such regions occur for r k ≤ 10 . The insignificant fluid motion at depth is consistent with our earlier studies based on the linear (Kameyama 2016) and nonlinear (Kameyama and Yamamoto 2018) analyses in Cartesian geometry.
The comparison of Fig. 3b and c with 3a also suggests that the changes in the values of ℓ c0 and K c0 with increasing r k for large r η in Fig. 2b can be understood from the changes in the vertical extent of incipient convection cells. For example, the increase in r k from 1 to 10 makes the stagnant lids along the top surface thinner and, hence, the convection cells beneath the lids taller. This in turn increases the horizontal extent of convection cells and, in other words, decreases ℓ c0 and K c0 of infinitesimal perturbations. Increasing r k from 10 to 100, by contrast, not only leads to thinner stagnant lids along the top surface but also thickens the regions of insignificant fluid motion above the bottom boundary. Owing to these combined effects, the vertical extent of convection cells becomes smaller as r k increases to 100, leading to an increase in ℓ c0 and K c0 of infinitesimal perturbations.
In Fig. 4, we show the variations in the radial profiles of eigenfunctions for radial velocity V for the critical states as functions of r η obtained for the various model geometries and several values of 1 ≤ r k ≤ 100 . Shown in colors and by thin contour lines in each plot of Fig. 4 are the values of V normalized by their maxima V max for given r η , while indicated by the green lines are the heights at which V = V max for given r η . In the cases where the regions of stable thermal stratification occur above the bottom boundary, we show by the thick gray lines the heights below which the thermal stratification is stable. Fig. 4 shows that a cold and stagnant lid develops along the top surface for sufficiently large r η , regardless of the values of r k and model geometries. Indeed, all the plots in Figure 4 show that the regions with very small V develop near the top surface when r η is sufficiently larger than the threshold around O (10 3 -10 4 ). In addition, for the cases with r k = 100 in 3D and 2D spherical geometries, the regions with very small V develop near the bottom boundary when r η is sufficiently large. The occurrence of such regions reflects the emergence of thick regions of stable thermal stratification above the bottom boundaries, and is consistent with the results of the Cartesian cases where such regions develop for r k ≥ 30. Moreover, from a detailed comparison of the location of the maxima in vertical velocity V between the cases with different r k in Fig. 4, we can see that the static stability of thermal stratification in the fluid layer exerts strong controls on the structures of incipient flows. For the cases where regions of stable thermal stratification are absent in the fluid layer ( r k ≤ 10 for 3D and 2D spherical geometries, and r k ≤ 3 for Cartesian geometry), the locations of the maxima in V shift downward with increasing r η , which approach to the bottom boundary for asymptotically large r η . For the other cases where regions of stable stratification exist, by contrast, the locations of the maxima in V, which shift downward with r η for small r η , do not reach the bottom boundary even when r η is asymptotically large. Instead, the asymptotic height of the Fig. 4 The radial profiles of the velocity perturbations V against r η in the fluid layer of different model geometries obtained for the cases with a r k = 1 , b r k = 3 , c r k = 10 , d r k = 30 and e r k = 100 . In each diagram, the vertical axis indicates the dimensionless height, while the horizontal axis indicates log 10 r η . Shown in color and thin solid contours are the amplitudes of radial velocity perturbations V normalized by their maxima V max for given r η . For the color scale for the distribution of V /V max , see the scale bar at the bottom of the figure. The contour lines are drawn for 0.1 ≤ |V /V max | ≤ 0.9 with the interval of 0.1. The green lines indicate the heights which yield V = V max for given r η . In some cases where regions of stable thermal stratification exist, we also show by thick gray lines the heights below which the thermal stratification is stable Kameyama Earth, Planets and Space (2021) 73:167 maxima in V for very large r η is equal to that at the top surface of the stable stratification indicated by the thick gray lines in the figure. In short, the incipient convection cells tends to shrink in the radial direction around the height of the top surface of the stable stratification with increasing r η , in response to the thickening of stagnant lids along the top surface and that of regions of insignificant fluid motion above the bottom boundary. We will conclude this section by demonstrating the similarity between the results of linear analysis and those of nonlinear one. In Fig. 5 we show the overall structure of stagnant-lid convection with large r k ( = 30 ) obtained in our earlier numerical experiments (Kameyama and Yamamoto 2018). Because of the large increase in thermal conductivity k with depth, the temperature gradient in the vertical direction is greatly reduced at depth (Fig. 5a,  b), and it eventually becomes smaller than the adiabatic gradient near the bottom boundary (Fig. 5c). This leads to an emergence of a region of stable thermal stratification at the bottom of the layer (Fig. 5d, e), although its dimensionless thickness ( ∼ 0.25 ) is smaller than that expected from the static stability of the reference conductive state ( ∼ 0.37 ; Fig. 1c). In this region, the fluid motion is indeed insignificant particularly in the vertical direction (Fig. 5f ). From this figure, we can confirm that the insignificant fluid motion in the basal region is largely due to the the stable thermal stratification there.

Conclusion
In this paper, we carried out a series of linear analyses on the onset of thermal convection of highly compressible fluids whose physical properties strongly vary in space in convecting vessels of either a three-dimensional (3D) spherical shell or a two-dimensional (2D) spherical annulus, by extending our earlier work in the Cartesian geometry (Kameyama 2016): the variations in thermodynamic properties (thermal expansivity and reference density) with depth are taken to be relevant for the super-Earths with ten times the Earth's mass, while the thermal conductivity and viscosity are assumed to exponentially depend on depth and temperature, respectively. One of the major finding in our analyses is that, for the cases with strong temperature dependence in viscosity and strong depth dependence in thermal conductivity, the critical Rayleigh number is on the order of 10 8 to 10 9 , which is close to the values of Rayleigh number of the thermal convection in the mantle of super-Earths with ten times the Earth's mass extrapolated from that in the Earth's mantle. This suggests that the mantle convection of massive super-Earths is most likely to fall in the stagnant-lid regime very close to the critical condition, if the properties of their mantle materials are quite similar to the Earth's.
Our analyses also demonstrated that the structures of incipient flows of stagnant-lid convection in the presence of strong adiabatic compression are strongly affected by the depth dependence in thermal conductivity and the geometries of convecting vessels, through the changes in the static stability of thermal stratification of the reference (conductive) state. For example, the increase in thermal conductivity with depth tends to enhance the stability at depth in the fluid layer, by raising the temperature and reducing the gradient in temperature in the radial direction at depth. When the increase is sufficiently large ( r k ≡ k bot /k top ≃ 100 ), regions of stable thermal stratification can occur in the deep mantles of super-Earths with a mass ten times that of the Earth. The occurrence of stable thermal stratification in the deep mantle strongly prevents the incipient cells of stagnant-lid convection from extending down to the bottom boundary, leading to the formation of "deep stratospheres" of insignificant fluid motions above the bottom hot boundaries in addition to the stagnant lids along the top cold surfaces. The occurrence of such motionless regions above the bottom boundaries is consistent with our earlier fully dynamic numerical experiments in Cartesian geometry (Kameyama and Yamamoto 2018). On the other hand, the spherical geometries of convecting vessels tend to suppress the occurrence of "deep stratospheres" above the bottom boundaries. This is because the reduction of surface area with depth in the spherical geometries steepen the gradient in temperature in the radial direction at depth. Furthermore, the increase in thermal conductivity k may not necessarily be so large as assumed here in silicate mantles of real super-Earths (e.g., Dalton et al. 2013;Hsieh et al. 2017Hsieh et al. , 2018Deschamps and Hsieh 2019), and the temperature dependence in k is likely to alter the nature of stagnant-lid convection (Deschamps 2021). Our analyses showed, however, that the regions of insignificant fluid motion can be formed above the bottom boundaries even in the three-dimensional spherical geometry in the presence of a factor of about 100 increase in the thermal conductivity with depth. We can, therefore, speculate that the stagnantlid convection in the mantles of massive super-Earths is accompanied by another motionless regions at the base of the mantles if the thermal conductivity of their mantle materials strongly increases with depth (or pressure).
The occurrence of stable thermal stratification in the mantles of massive super-Earths may strongly influence the mantle dynamics and surface environments of the particular terrestrial planets. For example, the stable region is most likely to greatly suppress the hot ascending plumes originating from the core-mantle boundary, resulting in the lack of hotspot volcanism on super-Earths. In addition, the stable region is most likely to slow down the cooling of underlying metallic cores, which may further imply the absence of intrinsic magnetic fields of super-Earths.
Our findings on the structures of incipient convection of highly compressible fluids may have strong impact on the studies on the evolution of massive terrestrial planets, in particular on the theoretical ones based on a parameterized convection (e.g. Tachinami et al. 2011;Foley 2019). In these earlier studies, the efficiency of convective heat transfer and emission of masses (such as greenhouse gases) from the planetary mantles are estimated, under the assumption that the convection structures in the mantles of massive super-Earths are similar to those of smaller terrestrial planets such as the Earth's. However, the occurrence of "deep stratospheres" presented in this study indicates that this assumption can be significantly violated in the mantles of massive terrestrial planets owing to the effect of adiabatic compression, in particular when the transport properties (i.e., viscosity and thermal conductivity) of mantle materials strongly vary in space. In addition, from the estimate of Rayleigh number and comparison with the critical value, the mantle convection of super-Earths is expected to have flow structures quite similar to the incipient ones presented here. We, therefore, conclude that the development of convection models which appropriately take into account the effect of strong adiabatic compression is of ultimate importance in deepening the understanding of the evolution of massive super-Earths.
where L 2 ≡ ℓ(ℓ + 1) . The equations for 2D spherical geometry are, on the other hand, In (23) and (25) we rewrote Ra with Ra c . The boundary conditions for are where r = r max and r = r min denote the top and bottom boundaries in nondimensional units, respectively. On the other hand, the boundary conditions for V are for 3D spherical geometry, and for 2D spherical geometry. We study the variation in the nature of the infinitesimal perturbations which become unstable for the smallest Rayleigh number by a following strategy. First, for given conditions (depth dependence in k and (24) (27) = 0 at r = r max , r min , temperature dependence in η ), we seek for a critical Rayleigh number Ra c of perturbations characterized by ℓ or K. In this procedure, we first recast the set of Eqs. (23) and (24) [or that of (25) and (26)] into an eigenequation of the form and then solved for an eigenvalue ( 1/Ra c ) together with eigenfunctions x (radial profiles of perturbations of temperature and radial velocity V) describing a critical state (i.e., growth rate a = 0 ). The numerical details on how to construct the matrix A can be found in Kameyama et al. (2015). The eigensystem developed thus was solved for all eigenvalues and eigenvectors by invoking the LAPACK subroutine "dgeev" designed for real nonsymmetric matrices. Among the series of eigenvalues and eigenfunctions in hand, we then search the ones which yield the absolute minimum of the critical Rayleigh number Ra c0 and the degree ℓ c0 or wavenumber K c0 corresponding to Ra c0 .
To conclude this appendix, we note that a special care must be taken in the final stage of the solution of the present eigensystem. Because the matrix A is in general a nonsymmetric matrix, the eigenequation (30) has to be solved in a field of complex numbers, even if the matrix A is real (e.g., Zebib et al. 1983;Zebib 1993). In fact, we obtained a set of complex eigenvalues and their conjugation. In addition, the distributions of eigenvalues of A significantly vary depending on whether a region of stable stratification (30) A · x = 1 Ra c x occurs or not. Indeed, all the eigenvalues have positive real parts (i.e. A is positive definite) when a thermal stratification is unstable in the entire layer. However, when a region of stable stratification occurs in the fluid layer, the eigenvalues have both positive and negative real parts (i.e. A is indefinite).  further showed that all the eigenvalues have negative real parts (i.e. A is negative definite) when a thermal stratification is stable in the entire layer.) Since the eigenvalues of A whose real parts are negative are physically unsound (see Kameyama et al. (2015) for detail), we can only focus on the eigenvalues 1/Ra c with positive real parts.