Linear Analysis On The Onset Of Thermal Convection of Highly Compressible Fluids With Variable Physical Properties In Spherical Geometry: Implications For The Mantle Convection of Super-Earths

In this paper we carried out a series of linear analysis on the onset of thermal convection of highly 7 compressible ﬂuids whose physical properties strongly vary in space in convecting vessels either of a 8 three-dimensional spherical shell or a two-dimensional spherical annulus. The variations in 9 thermodynamic properties (thermal expansivity and reference density) with depth are taken to be 10 relevant for the super-Earths with 10 times the Earth’s mass, while the thermal conductivity and 11 viscosity are assumed to exponentially depend on depth and temperature, respectively. Our analysis 12 showed that, for the cases with strong temperature-dependence in viscosity and strong 13 depth-dependence in thermal conductivity, the critical Rayleigh number is on the order of 10 8 to 10 9 , 14 implying that the mantle convection of massive super-Earths is most likely to fall in the stagnant-lid 15 regime very close to the critical condition, if the properties of their mantle materials are quite similar to 16 the Earth’s. Our analysis also demonstrated that the structures of incipient ﬂows of stagnant-lid 17 convection in the presence of strong adiabatic compression are signiﬁcantly aﬀected by the 18 depth-dependence in thermal conductivity and the geometries of convecting vessels, through the changes 19 in the static stability of thermal stratiﬁcation of the reference state. When the increase in thermal 20 conductivity with depth is suﬃciently large, the thermal stratiﬁcation can be greatly stabilised at depth, 21 further inducing regions of insigniﬁcant ﬂuid motions above the bottom hot boundaries in addition to 22 the stagnant lids along the top cold surfaces. We can therefore speculate that the stagnant-lid 23 convection in the mantles of massive super-Earths is accompanied by another motionless regions at the 24 base of the mantles if the thermal conductivity strongly increases with depth (or pressure), even though 25 their occurrence is hindered by the eﬀects the spherical geometries of convecting vessels. 26

Super-Earths are extra-solar planets which have small masses (compared to that of Jupiter) and high 31 mean density larger than 5000 kg/m 3 (Howard et al. 2010). Since the first discovery of such massive 32 terrestrial planets (Rivera et al. 2005), a growing effort has been devoted to the studies in the internal 33 (radial) structure (e.g., Seager et al. 2007;Valencia et al. 2007b;Sotin et al. 2007;Wagner et al. 2011), 34 the thermal evolution (e.g., Papuc and Davies 2008;Kite et al. 2009;Korenaga 2010;Gaidos et al. 2010;35 Tachinami et al. 2011;Čížková et al. 2017) and the habitability (e.g., Noack et al. 2017;Tosi et al. 2017;36 Foley 2019) of super-Earths. During the last decades, the studies on mantle dynamics have been playing 37 an important role in addressing these issues. 38 It can be intuitively understood that, because of their large sizes, the convection in the mantles of massive 39 super-Earths is most likely to occur in a quite different manner from those in the Earth and smaller 40 terrestrial planets. For example, the variations of physical properties of mantle materials are expected to 41 be larger in the much thicker mantles of super-Earths than in the Earth's (e.g., Karato 2011). Indeed, 42 the effects of strong spatial variations in physical properties (such as viscosity and thermal conductivity) 43 on the convecting flow patterns have been widely studied so far, under the assumption that the effects 44 of other physical mechanisms such as compressibility can be minor (e.g., Kameyama and Ogawa 2000;45 Miyauchi and Kameyama 2013). In particular, the effects of strong temperature-dependence in viscosity 46 have been intensively studied by numerical models of convection in three-dimensional spherical geometry 47 (e.g., Yanagisawa et al. 2016;Guerrero et al. 2018). 48 On the other hand, it has been well acknowledged that the effect of adiabatic compression affects the 49 nature of thermal convection (Jarvis and McKenzie 1980;King et al. 2010;Liu and Zhong 2013). This  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 in physical properties of mantle materials on the convecting flow structures in spherical vessels, in order 60 to gain the fundamental insights into the mantle convection of massive super-Earths. 61 In this study, we carry out a series of linear analysis on the thermal convection in spherical domains of a 62 highly compressible fluid with the spatial variation in physical properties, in order to investigate how the 63 spherical geometries of planetary mantles can affect the critical state and condition of thermal convection 64 in the presence of the interplay between the adiabatic compression and spatial variations in physical 65 properties such as viscosity and thermal conductivity. In particular, we will focus on the influences of The onset of thermal convection is considered for a highly compressible fluid of infinite Prandtl number, 72 as a model of mantle convection of a super-Earth with 10 times the Earth's mass. The numerical model 73 is the same as those of our earlier one (Kameyama 2016) except that the convecting vessels are taken 74 to be either a three-dimensional (3D) spherical shell or a two-dimensional (2D) spherical annulus whose 75 inner and outer radii are r min and r max , respectively, and r min /r max = 0.5. The temperature T of the 76 fluid layer is fixed to be T = T top and T bot (> T top ) at the top and bottom boundaries, respectively. In 77 order to keep the present analysis as simple as possible, we ignored the effects of internal heating and 78 spatial variations in specific heat C p of the fluid and the gravity g throughout the layer. In addition we 79 fixed T top /∆T = 0.1 (where ∆T ≡ T bot − T top is the temperature difference across the layer), indicating 80 that the temperature at the top surface of the planet is similar to the Earth's.

81
Physical properties of modelled compressible fluid 82 As a model of pressure-dependence in thermodynamic properties of mantle materials of massive super-Earths, we employ the same variations in thermal expansivity α and reference density ρ with depths as in our earlier works (Tachinami et al. 2014;Miyagoshi et al. 2014Miyagoshi et al. , 2015Kameyama 2016;Kameyama and Yamamoto 2018). The thermal expansivity α decreases with depth, whose 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 10 times the Earth's mass. That is, the value of α(r = r min ) at the bottom 4   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 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 83 of ρ(r = r min ) at the bottom boundary larger than ρ(r = r max ) ≡ ρ top at the top surface by about 84 3.17 (Valencia et al. 2006). We note that the linear increase in ρ assumed here is not likely to spoil 85 the significance of the present analysis despite its geophysical unsoundness, since the vertical profile of 86 (reference) density has a very minor effect on the nature of critical states of thermal convection of highly 87 compressible fluids (Kameyama et al. 2015).

88
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 given by 89 r k = k(r = r min )/k(r = r max ). For a geophysical significance, we restrict ourselves to the cases where 90 r k ≥ 1; the thermal conductivity k increases with depth for r k > 1. According to the recent estimates super-Earths. In this study we varied r k up to 100.

95
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 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, of ∆T ≡ T bot − T top . Here ⊗ stands for the tensor product of vectors, v ⊗ ∇ is the transpose of 108 ∇ ⊗ v, I is the identity tensor, e r is the unit vector in radial direction positive outward, v is the 109 velocity vector, v r ≡ (v · e r ) is the velocity in the radial direction, p is pressure, andε II is the second 110 invariant of strain rate tensor. In the right-hand side of (6), the first, second, and third terms represent 111 the effects of apparent heat transport (advection plus conduction), adiabatic temperature change, and 112 viscous dissipation (frictional heating), respectively.

113
The present formulation includes two nondimensional parameters. First one is the dissipation number Di of the modelled fluid defined by, Second one is the Rayleigh number Ra of thermal convection defined by, of super-Earths whose mass is 10 times larger than the Earth's (Tachinami et al. 2014;Miyagoshi et al. 121 2014Miyagoshi et al. 121 , 2015Miyagoshi et al. 121 , 2018. In the present study, we fixed Di = 5, while the values of Ra estimated above can be 122 compared with those of the critical Rayleigh numbers derived from the linear stability analysis.

123
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 124 is the reason why the effect of viscous dissipation is eliminated in (14).

125
At the top and bottom boundaries, we kept the perturbations in temperature T ′ and the radial velocity 126 v ′ r to be zero. For the horizontal velocity at the boundaries, on the other hand, we assumed a free-slip 127 condition.

128
Numerical method

129
The equations for infinitesimal perturbations (12) to (14) are recast into the one-dimensional differential 130 equations in the radial direction, by applying the spectral expansion in the horizontal planes using 131 either the spherical harmonics or sinusoidal function depending on three-and two-dimensional spherical 132 geometries. The derived one-dimensional differential equations are discretised with an eighth-order 133 accuracy in space using equally-spaced 257 grid points. The eigenvalue problem for a given condition 134   7   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 is more significant in 3D cases (∝ r −2 ) than in 2D cases (∝ r −1 ). Owing to this effect, the magnitude of 161 |dT /dr| in the shallow part of fluid layer is smaller in 3D cases than in 2D cases and, at the same time, 162 the magnitude in the deep part is larger (see Figure 1b). This results in lower T in the entire layer for 3D presented in Figure 1, indeed, the values of T for given r k are the highest for Cartesian cases where the 166 effect of the changes in the surface area with depth is absent.

167
In order to see the static stability of the thermal stratification more clearly, we show in Figure 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 grey 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 (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 168 static stability of thermal stratification is stronger for larger N 2 and dT pot /dr.

Availability of data and materials 322
The data that support the findings of this study and the numerical programs used in this study the are 323 available upon reasonable request to M.K. (email: kameyama@sci.ehime-u.ac.jp).

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

Figure 3
Schematic illustrations on two-dimensional planes of the structures of incipient ows obtained for the cases with strong temperature-dependent viscosity rη = 108 and various model geometries together with the values of (a) rk = 1, (b) rk = 10 and (c) rk = 100. Note that for the cases with 3D spherical geometry we showed the ow structures in meridian planes assuming the incipient ows occurring with only zonal components (i.e., the spherical harmonic order m = 0). Shown in colour are the distributions of in nitesimal perturbations in temperature T0 normalised by their maxima Θmax. For the colour scale for the distribution of T=Θmax, see the scale bar at the bottom of the gure. Shown in solid lines are the contour lines of a "potential" Ψ of the mass ux pv'. The contour lines are drawn for 0:1 ≤ | Ψ/Ψmax| ≤ 0:9 with the interval of 0.1, where max is the maximum of |Ψ|.

Figure 4
The diagrams for the radial pro les of the velocity perturbations V against rη in the uid layer of different model geometries obtained for the cases with (a) rk = 1, (b) rk = 3, (c) rk = 10, (d) rk = 30 and (e) rk = 100. In each diagram, the vertical axis indicates the dimensionless height, while the horizontal axis indicates log10 rη. Shown in colour and thin solid contours are the amplitudes of radial velocity perturbations V normalised by their maxima Vmax for given rη. For the colour scale for the distribution of V=Vmax, see the scale bar at the bottom of the gure. The contour lines are drawn for 0.1 ≤ |V/Vmax| ≤ 0.9 with the interval of 0.1. The green lines indicate the heights which yield V = Vmax for given rη. In some cases where regions of stable thermal strati cation exist, we also show by thick grey lines the heights below which the thermal strati cation is stable.