# Compressible convection in super-earths

Yanick Ricard, Thierry Alboussière

► To cite this version:

Yanick Ricard, Thierry Alboussière. Compressible convection in super-earths. Physics of the Earth and Planetary Interiors, 2023, 341, pp.107062. 10.1016/j.pepi.2023.107062. hal-04289858

HAL Id: hal-04289858

https://hal.science/hal-04289858v1

Submitted on 16 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Compressible convection in super-earths

Yanick Ricard*, Thierry Alboussière

Université de Lyon, ENSL, UCBL, UJM, Laboratoire LGLTPE, 15 parvis René Descartes, BP7000, 69342, Lyon Cedex 07, France

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

**Keywords:**

Super-Earth
Compressible convection
Marginal stability
Adiabatic conditions

## ABSTRACT

The radial density of planets increases with depth due to compressibility, leading to impacts on their convective dynamics. To account for these effects, including the presence of a quasi-adiabatic temperature profile and entropy sources due to dissipation, the compressibility is expressed through a dissipation number, $\mathcal{D}$, proportional to the planet's radius and gravity. In Earth's mantle, compressibility effects are moderate, but in large rocky or liquid exoplanets (super-earths), the dissipation number can become very large. This paper explores the properties of compressible convection when the dissipation number is significant. We start by selecting a simple Murnaghan equation of state that embodies the fundamental properties of condensed matter at planetary conditions. Next, we analyze the characteristics of adiabatic profiles and demonstrate that the ratio between the bottom and top adiabatic temperatures is relatively small and probably less than 2. We examine the marginal stability of compressible mantles and reveal that they can undergo convection with either positive or negative superadiabatic Rayleigh numbers. Lastly, we delve into simulations of convection in 2D Cartesian geometry performed using the exact equations of mechanics, neglecting inertia (infinite Prandtl number case), and examine their consequences for super-earth dynamics.

## 1. Introduction

The interior of the Earth's mantle and that of other solid planets of the solar system are driven by convection. These motions are studied numerically using convection codes that often ignore the compressibility of the material and use the Boussinesq approximation in which the density of the fluid is uniform except for the buoyancy term (see e.g., Choblet et al., 2007; Zhong et al., 2008). In some cases, anelastic equations are used that consider first-order perturbations around an adiabatic state that accounts for the change of parameters with depth (see e.g., Zhang and Yuen, 1996; Tackley, 2008; King et al., 2010; Kameyama and Yamamoto, 2018). These anelastic codes often used a simplified equation of state (the relationship between density, temperature, and pressure) and implement some approximations that contradict the thermodynamic rules: this happens most often when some thermodynamic quantities are given an explicit dependence in depth rather than in pressure and temperature. We have shown in a series of papers that, when inertia can be neglected, which is the case in the mantle of planets, solving the exact mechanical and thermodynamical equations for a realistic equation is however not much more difficult than implementing an anelastic formalism (Alboussière and Ricard, 2017; Alboussière et al., 2022; Curbelo et al., 2019; Ricard et al., 2022).

Large exoplanets have been detected with densities that suggest that their composition is probably similar to the telluric planets in our system. While compressibility effects remain modest in the Earth's mantle and core, they become an essential ingredient of Super-Earths' internal dynamics. In this paper we study the Rayleigh-Bénard convection when the dissipation number is large, by solving the exact set of compressible equations without inertia. In a first section we choose an appropriate equation of state for condensed materials, in a second section we discuss some generic characteristics of the adiabatic conditions, in a third section we discuss the values of the critical Rayleigh number at the onset of convection as a function of the dissipation number and the imposed temperature difference and we present various 2D Cartesian cases of developed compressible convection depending on these two parameters.

## 2. Adiabatic conditions inside a convective planet

It is well known (see e.g., Braginsky and Roberts, 1995; Schubert et al., 2001; Ricard, 2015) that convection of a compressible fluid at high Rayleigh number brings the average radial profiles of density, temperature and pressure close to their adiabatic and hydrostatic values ($\rho_a$, $Ta$, $P_a$) according to

<footnote>* Corresponding author.
E-mail address: ricard@ens-lyon.fr (Y. Ricard).</footnote>
$$ \frac{d\ln\rho_a}{dz} + \frac{\alpha_a g}{\Gamma_a C_p^a} = 0, \qquad (1a) $$

$$ \frac{d\ln T_a}{dz} + \frac{\alpha_a g}{C_p^a} = 0, \qquad (1b) $$

$$ \frac{dP_a}{dz} + \rho_a g = 0, \qquad (1c) $$

z being the vertical coordinate (directed against gravity g = -gez). In these equations, α is the thermal expansivity, Cₚ is the heat capacity (or specific heat) at constant pressure and Γ is the Grüneisen parameter

$$ \Gamma = \frac{\alpha K_T}{\rho C_V} = \frac{1}{\rho C_V} \left(\frac{\partial P}{\partial T}\right)_V. \qquad (2) $$

The superscript or subscript 'a' in eqs. (1a)-(1b)-(1c) indicates that the various quantities are computed along the adiabat itself.

A reasonable equation of state (EoS) for a condensed planet is based on the observation that the Grüneisen parameter is essentially a function of density (Anderson, 1979) according to

$$ \Gamma = \Gamma_0 \left(\frac{\rho_0}{\rho}\right)^q, \qquad (3) $$

where q is around 1, ρ₀ and Γ₀ are the density and the Grüneisen parameter at standard conditions that we choose to be on the surface of the planet. On average, the Grüneisen parameter is between 1 and 2 in the mantle (e.g. Stacey and Davis, 2004) or in the core (Alfe et al., 2002). In what follows, we will use q = 1. With the additional assumption that C_V is constant, the EoS is necessarily of the form P = f(ρ) + α₀K_T^0(T - T₀), i.e. a function of density plus a thermal contribution where α₀ and K_T^0 are the thermal expansivity and isothermal incompressibility under reference conditions. For condensed materials, a suitable function for the pressure-density relationship at constant temperature is given by a Murnaghan expression (Murnaghan, 1951) (a special case of a polytropic EoS, (Chandrasekhar, 1957)),

$$ P = \frac{K_T^0}{n} \left[ \left( \frac{\rho}{\rho_0} \right)^n - 1 \right] + \alpha_0 K_T^0 (T - T_0), \qquad (4) $$

where n ≈ 3–4 for solid silicates and for liquid silicates or metals. Although EoS (4) is simple and empirical, it encompasses the typical solid and fluid properties and gives a very good fit to the Earth's radial density assuming its adiabaticity, away from the transition zone discontinuities (Ricard et al., 2022). This EoS implies direct relationships between thermal expansivity and incompressibility with density that are

$$ \alpha = \alpha_0 \left(\frac{\rho_0}{\rho}\right)^n, \qquad (5) $$

$$ K_T = K_T^0 \left(\frac{\rho}{\rho_0}\right)^n, \qquad (6) $$

and again these two expressions provide realistic expressions of the properties measured in laboratory experiments. Since the Grüneisen parameter (3) depends only on density, we obtain a simple relation between the adiabatic temperature and the adiabatic density by combining (1a) and (1b)

$$ T_a = T_a^t \exp \left[ \Gamma_0 \left( \frac{\rho_0}{\rho_a^t} - \frac{\rho_0}{\rho_a} \right) \right], \qquad (7) $$

where T_a^t and ρ_a^t are the surface (t stands for top) adiabatic temperature and density.

We make two additional approximations when deriving analytical expressions for the adiabatic profiles (those approximations will not be used in the numerical computations solving the exact equations).

* We assume the heat capacities at constant volume and temperature C_V and C_P are equal and constant. The heat capacites are indeed very close when αT << 1 which is the case for condensed materials, and they are both close to the Dulong and Petit values C_V ≈ C_P ≈ 3R in J K⁻¹ mol⁻¹ (R is the gas constant) at large temperatures (Petit and Dulong, 1819).

* The surface adiabatic temperature is different from the surface temperature, but since α(T_a^t - T₀) < 1, the surface adiabatic density and the reference density can also be identified, ρ_a^t ≈ ρ₀ (in other words the density of planets is mainly a function of pressure, not temperature).

A last assumption is made on the variation of gravity with depth in a generic planet. For simplicity, we assume that gravity is uniform which is basically the case in Earth's mantle.

With these hypotheses it is easy to solve for the adiabatic conditions in a layer where z varies between 0 and H (e.g., a mantle of thickness H, z = 0 being at the core-mantle boundary, see also Zhang and Yuen (1996) starting with a slightly different EoS), and we get

$$ \rho_a = \rho_0 \left( 1 + \frac{H-z}{h} \right)^{1/(n-1)}, \qquad (8a) $$

$$ P_a = \frac{n-1}{n} \rho_0 g h \left[ \left( \frac{\rho_a}{\rho_0} \right)^n - 1 \right], \qquad (8b) $$

$$ T_a = T_a^t \exp \left[ \Gamma_0 \left( 1 - \frac{\rho_0}{\rho_a} \right) \right] \qquad (8c) $$

where

$$ h = \frac{1}{n-1} \frac{K_T^0}{\rho_0 g} = \frac{1}{n-1} \frac{\Gamma_0}{\mathcal{D}} H. \qquad (9) $$

In the last equality, we have introduced the dissipation number D defined by

$$ \mathcal{D} = \frac{\alpha_0 g H}{C_V}. \qquad (10) $$

This surface dissipation number is only expressed from quantities known at the surface, this choice seems to be the only possible choice when exploring a new planet. In the Earth the dissipation number is around D⊕ = 0.71 in the mantle and 0.56 in the liquid core (using α₀ = 3 × 10⁻⁵ K⁻¹, H = 2900 km and C_V = 1200 J K⁻¹ kg⁻¹ in the mantle, α₀ = 1.8 × 10⁻⁵ K⁻¹ (Murphy et al., 2013), H = 2300 km and C_V = 715 J K⁻¹ kg⁻¹ (Gubbins et al., 2003) in the liquid core, with g = 9.8 m s⁻²).

In geophysical textbooks (see, e.g., Schubert et al., 2001) D is defined with a Cₚ in the denominator which does not make much practical difference as their difference is always neglected in the geology literature. However we prefer to define the dissipation number with C_V as in the definition of the Grüneisen parameter. In order to obtain a simple expression for the EoS, we have made the choice of a constant C_V, hence Cₚ cannot be a constant too because their difference has to obey Mayer's relation. Assuming the constancy of both heat capacities would lead to inconsistencies in the energy conservation (Alboussière and Ricard, 2013).

The preceding Eqs. (8a)-(8b)-(8c) can be used to discuss possible characteristics of the adiabatic profiles of large planets. From the variety of masses and radii of exoplanets that have been detected, it appears that many of them are rocky at least until a radius of order 2.5 times the Earth's radius (Otegi et al., 2020). Their observed mass M increases roughly as a power 3.45 of their radius R (their large internal pressures increase their average densities as ≈ R⁰.⁴⁵). We will use this observation to scale the gravity in our equations with gM/R² ≈ R¹.⁵ and consider that the thickness of the convective layers is proportional to R. With these scalings, the dissipation number D varies like gH/R².⁵. Therefore,
according to Otegi et al. (2020) dissipation numbers up to 2.5<sup>2.5</sup> = 10 times that of the Earth can be expected in fairly common rocky planets, and in what follows we will explore dissipation numbers up to ℳ = 10. We will use ℳ = ℳ<sub>⊕</sub>(R/R<sub>⊕</sub>)<sup>2.5</sup> when, to fix ideas, we discuss in term of planetary radii instead of dissipation numbers; a Super-Earth with a radius twice that of the Earth (resp. 3 times) would therefore be assumed to have a mantle with a dissipation number around 4.0 (resp. 11.1) and a core with a dissipation number around 3.2 (resp. 8.7).

### 3. Adiabatic conditions in a Super-Earth

#### 3.1. The adiabatic density and temperature profiles

Since incompressibility increases very significantly with density and thus with pressure, the adiabatic density and temperature increase only moderately with planetary radius. In Fig. 1, we plot the adiabatic temperature normalized with its value at the surface according to (8c). We use ℳ = ℳ<sub>⊕</sub> = 0.6 (black), ℳ = 2 (red), ℳ = 10 (green).

The ratio of the adiabatic density between the bottom and the top of a convecting mantle is according to (8a) and (9)

$$ \frac{\rho_a^b}{\rho_a^t} = \left(1 + (n-1) \frac{\mathcal{D}}{\Gamma_0}\right)^{\frac{1}{n-1}}. \qquad (11) $$

This ratio is plotted in Fig. 2a as a function of ℳ (bottom axis) and as a function of R (top axis, assuming ℳ ∝ R<sup>2.5</sup>). The ⊕ symbol indicates the position of the Earth where an adiabatic density ratio of 1.52 is predicted through the mantle (due to the phase changes in the transition zone the observed density change in Earth's mantle is rather 1.70).

This adiabatic density ratio controls the adiabatic temperature ratio according to (8c) (see Fig. 2b). For the Earth, this ratio should be 1.41 through the mantle (indicated by the symbol ⊕, say from 1600 K on top to 2256 K at the bottom). The maximum bottom temperature T<sub>a</sub><sup>b</sup> is, at any rate, bounded when ρ<sub>a</sub> → ∞ by

$$ \frac{T_a^b}{T_a^t} \le e^{\Gamma_0} \approx 2.72. \qquad (12) $$

Even in very large silicated Super Earth, the bottom adiabatic temperature should remain moderate and should hardly exceed a factor of 2 times the surface adiabatic temperature (Fig. 2b).

<figure><img src="image_3.png" /><figcaption>Fig. 1. Normalized adiabatic temperature for ℳ = 0.6, 2, 10 (black, red, green). In the case ℳ = 10, the dashed lines are conductive profiles that will be discussed below, one is with ΔT = ΔT<sub>a</sub> (dark green), one with a gradient equal to the bottom adiabatic gradient (orange), the last one (blue) carries at the surface the same heat flow as the adiabatic gradient at the surface. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

#### 3.2. The adiabatic temperature gradient

According to (1b), the surface adiabatic gradient is simply

$$ \frac{dT_a}{d\tilde{z}}|_t = \mathcal{S} T_a^t, \qquad (13) $$

where $\tilde{z} = z/H$ is the normalized height in the convective layer. Obviously the heat carried out near the surface, along the adiabat, increases with ℳ. However the adiabatic gradient near the bottom is

$$ \frac{dT_a}{d\tilde{z}}|_b = \mathcal{S} \left( \frac{\rho_a^t}{\rho_a^b} \right)^n T_a^b. \qquad (14) $$

In this expression, T<sub>a</sub><sup>b</sup> is bounded by eq. (12) and the thermal expansivity (related to the (ρ<sub>a</sub><sup>t</sup>/ρ<sub>a</sub><sup>b</sup>)<sup>n</sup> term) decreases faster than 1/S. This means that the adiabatic gradient at depth (in absolute value), initially increases and then decreases with the dissipation number (inspection of (11) shows that the adiabatic gradient at depth decreases with ℳ<sup>-1/(n-1)</sup> ≈ ℳ<sup>-0.43</sup>). This is visible in Fig. 1 (compare the adiabatic gradient at the bottom of the three curves). This implies that, as the dissipation number increases, the effects of compressibility become confined to shallow depths, while at depth, the fluid appears more and more incompressible. Unexpectedly, as the effects of compression increase (as the planet radius increases), deep convection appears more and more incompressible!

Another way to understand this is to consider that the compressible effects that affect convection are not related to ℳ, which is based on the reference thermal expansivity but to ℳ = ∫<sub>0</sub><sup>H</sup> ℳd\tilde{z}, where the dissipation number (10) is averaged over the layer thickness. Using the expressions (5) and (8a), as shown in Ricard et al. (2022), one gets

$$ \mathcal{S} \le \Gamma_0, \qquad (15) $$

i.e., the average dissipation number is never larger than the Grüneisen parameter which is about 1. The weak variation of the average dissipation number with planetary size was also noticed by Miyagoshi et al. (2015).

### 4. Convection in the mantle of super-Earths

#### 4.1. Compressible convection

In situations where compressibility is important and the physical parameters vary strongly with depth, the use of a simple Boussinesq convection model and correction of the results by adding an adiabatic contribution a posteriori are not sufficient. Using anelastic formulations (Ogura and Phillips, 1962; Jarvis and McKenzie, 1980; Braginsky and Roberts, 1995; Lantz and Fan, 1999) may also be tricky, as it is easy to inadvertently contradict the basic thermodynamic rules (see e.g., Leng and Zhong, 2008; Alboussière and Ricard, 2013). In a previous paper (Ricard et al., 2022), we explained how to solve the fully compressible equations without approximations, when inertia is neglected (the infinite Prandtl number approximation), which is appropriate for mantle convection, i.e., how to solve the following equations

$$ \frac{D\rho}{Dt} + \rho \nabla \cdot \mathbf{u} = 0, \qquad (16a) $$

$$ \eta \nabla^2 \mathbf{u} + \frac{\eta}{3} \nabla \cdot \nabla \mathbf{u} - \nabla P + \rho \mathbf{g} = 0, \qquad (16b) $$

$$ \rho T \frac{D \mathcal{S}}{Dt} = \dot{\mathbf{e}} : \mathbf{\tau} + k \nabla^2 T, \qquad (16c) $$

where the viscosity η and thermal conductivity k are assumed to be uniform. Following exactly the rules of thermodynamics and starting from the EoS (4), the entropy can be expressed and by integration of TdS = C<sub>V</sub>dT - αK<sub>T</sub>Tdρ/ρ<sup>2</sup> and writes
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Ratio between bottom and top adiabatic density (panel a) and temperature (panel b) in a compressible planet as a function of dissipation number $\mathcal{S}$ ($n = 3.3$, $\Gamma_0 = 1$). The symbol $\oplus$ indicates the situation for the Earth. The horizontal axis is either labelled in dissipation numbers (bottom) or in Super-Earth radii (top).</figcaption></figure>

$$ \mathcal{S} = C_V \ln \frac{T}{T_a} + \alpha_0 K_T^0 \left( \frac{1}{\rho} - \frac{1}{\rho_a} \right). \qquad (17) $$

which cancels, at it should, when the density and temperature are those of the adiabatic conditions. The adiabatic density and temperature are calculated from (1a) and (1b) (see Ricard et al., 2022, for details) where the heat capacity at constant pressure that appears in the adiabatic profile is exactly given by Mayer's relation

$$ C_P = C_V \left[ 1 + \Gamma_0 \alpha_0 T \left( \frac{\rho_0}{\rho} \right)^{n+1} \right]. \qquad (18) $$

## 4.2. Marginal stability and Schwarzschild criterion

When heated from below, a fluid begins to convect when two conditions are realized. First the local temperature gradient $|dT/dz|$ must be larger than the adiabatic gradient $|dT_a/dz|$. This is the Schwarzschild criterion (Schwarzschild, 1906): the temperature of a rapidly moving upward parcel of fluid follows the adiabatic gradient and must become warmer (i.e., less dense) than the surrounding medium to be gravitationally unstable. This criterion defines a necessary condition for convection. Second, the total temperature drop $\Delta T = T^b - T^c$ across the convective layer must be large enough so that a dimensionless number, the Rayleigh number, exceeds some critical value. In the simple case where the top and bottom boundaries are free slip, Rayleigh (1916) proved that this last condition can be expressed as

$$ Ra^c = \frac{\alpha_0 \rho_0^2 C_P g H^3 \Delta T}{\eta k} \geq \frac{27}{4} \pi^4 = 657.51 \qquad (19) $$

This sufficient condition for convection was obtained in the Boussinesq approximation where all the parameters $\alpha_0$, $\rho_0$, $\eta$ and $k$ are uniform. In a compressible fluid $\alpha$ and $\rho$ are depth dependent and (19) is not necessarily in agreement with the Schwarzschild criterion. As we are accustomed to thinking of the adiabatic gradient as more or less uniform (this would be exactly true if the fluid were a perfect gas and which is approximately true in the Earth's mantle), a superadiabatic Rayleigh number $Ra_{sa}$ is usually defined where the temperature drop $\Delta T$ is replaced by the superadiabatic temperature drop $\Delta T_{sa} = \Delta T - \Delta T_a$, the temperature drop in excess of the adiabatic temperature drop

$$ Ra_{sa} = \frac{\alpha_0 \rho_0^2 C_P g H^3 \Delta T_{sa}}{\eta \kappa} = Ra \frac{\Delta T_{sa}}{\Delta T}, \qquad (20) $$

where $\alpha_0$ and $\rho_0$ are now some characteristic values of the depth-dependent thermal expansivity and density. With this definition, $Ra_{sa}^c \geq$

657.51 is in agreement with both what is found in the Boussinesq approximation and with the Schwarzschild criterion (Malkus, 1954; Jeffreys, 1930; Grossmann and Lohse, 2001), at least, if $T_a$ is assumed to vary linearly with depth (i.e., if $dT_a/dz$ is uniform with $dT_a/dz = -\Delta T_a/H$). Various relations obtained in the Boussinesq approximation, for example between heat flow and Rayleigh number, are often taken to be valid in the compressible case when the superadiabatic Rayleigh number is used. We will show that the situation is more complex in super-Earth cases where the adiabatic gradient is large and with a large curvature.

To begin with a simple case where the adiabatic gradient is uniform, we consider a perfect gas with EoS $P = \rho RT$. This EoS is certainly not appropriate for a planetary mantle but it corresponds to the prototypical case of compressible convection. For a perfect gas, $\alpha T = 1$ and $C_P$ is constant in (1b) so that $dT_a/dz = -\alpha_a T_a g / C_P = -g / C_P$. Using the surface temperature $T_0$ and the height $H$ of the convection layer to non-dimensionalize the variables, we obtain $dT_a/dz = -\mathcal{S}/\gamma$ where $\gamma = C_P / C_V$ is the heat capacity ratio (also known as adiabatic index or Laplace's coefficient). A bottom temperature $rT_0$ is then considered to be imposed, in which case the conductive geotherm is $dT_c/dz = -(r-1)$. The Schwarzschild criterion imposes therefore that convection cannot exist when $\mathcal{S} \geq \gamma(r-1)$.

In a previous paper, we give the general equations verified by the marginally stable solution and how to compute the critical Rayleigh number for any EoS (eqs. 5.5–5.8 in Alboussière and Ricard (2017) and following comments). We therefore calculate the critical Rayleigh number $Ra^c(\mathcal{S}, r)$ for Rayleigh Bénard convection of a perfect gas and plot the result in Fig. 3a. As expected, convection can only occur below the $\mathcal{S} = \gamma(r-1)$ line. The cyan line corresponds to $Ra^c = 657.51$, and indeed for $\mathcal{S} \to 0$ and $r \to 0$, the critical value obtained for the Boussinesq case is recovered. Increasing the temperature jump $r-1$ decreases the critical Rayleigh number, increasing the dissipation number $\mathcal{S}$ increases the critical Rayleigh number.

The situation is obviously very different in a planet with a Murnaghan EoS. The motionless diffusive and hydrostatic state still has a linear geotherm $dT_c/dz = -(r-1)$ but the adiabatic gradient appearing in the Schwarzschild criterion depends on depth (see also Kameyama and Yamamoto, 2018). This adiabatic gradient is $dT_a/dz = -\mathcal{S} \rho_c^{-n} T_c / \gamma_c$, i.e., the slope of the adiabat under the conditions corresponding to the diffusive state with temperature $T_c$, density $\rho_c$ and heat capacity ratio $\gamma_c$. We have tested the conductive profiles over an extensive range of input parameters $\alpha_0 T_0$, $\Gamma_0$, $r$, $\mathcal{S}$, $n$, and observed that the minimum value of the adiabatic gradient is attained at the bottom of the layer in most cases. In other cases, the minimum is at the top of the layer, only when the dissipation number is small and the temperature ratio large, so that
<figure><img src="image_5.png" /><figcaption>Fig. 3. Isocontours of the critical superadiabatic Rayleigh number as a function of the surface dissipation number $\mathcal{S}$, and the ratio between top and bottom temperatures $r$. In the left panel, the EoS of the convective fluid is that of a perfect gas (with $\gamma = 5/3$), in the right panel we consider the Murnaghan EoS, appropriated for condensed matter. While for the perfect gas, a large dissipation number $\mathcal{S}$ decreases the domain where convection can occur, for the Murnaghan solid, convection can easily start in the deep layers even when the imposed temperature difference is lower than the adiabatic temperature difference. The dashed green and orange lines corresponds to the conductive slopes of Fig. 1. The cyan lines in both panel are for $Ra_{sa}^c = 657.51$, the value obtained by Lord Rayleigh in the Boussinesq case. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the adiabatic gradient is smaller than the conduction gradient all along the vertical profile (this can also be analytically proven if one assumes that the adiabatic gradient cannot be much different from $dT_a/dz \approx -\mathcal{S}\rho_a^n T_c$, i.e., that $\rho_c \approx \rho_a$ (see (8a)) and $\gamma_c \approx 1$). This implies that convection cannot occur when the diffusive gradient $|dT/dz|$ is lower than that of the dashed orange line in Fig. 1 (tangent to the adiabatic profile at the bottom where the adiabatic gradient is minimal in absolute value). When the adiabatic gradient is smaller than the conduction gradient only in the lower part, convection can start in the deep layers although $\Delta T \le \Delta T_a$, i.e., although the superadiabatic Rayleigh number is negative. This is confirmed by the computation of the critical Rayleigh number shown in Fig. 3b. There is a large domain in blue, where the conductive temperature gradient is in between that of the dark green curve and that of the orange curve of Fig. 1, where convection can start in the deep layer while the critical superadiabatic Rayleigh number is negative. Notice that again, when $\mathcal{S} \to 0$ and $r \to 0$, the critical value computed for the Boussinesq case is recovered ($Ra_{sa}^c \to 657.51$, a value shown by a cyan line). This limit is indeed independent of the chosen EoS.

A negative superadiabatic Rayleigh number is surprising at first sight. We could have discussed our results in terms of the total Rayleigh number only (see (20)), which remains positive or define another superadiabatic Rayleigh number that will necessarily be positive when the fluid convects (e.g., we could replace $\Delta T_a$ in (20) by $-H(dT_a/dz)(z=0)$, so that a convective Schwarzschild criterion is reached at the bottom only when the superadiabatic Rayleigh number is positive). However, changing the expression of the superadiabatic Rayleigh number each time a new EoS is considered, seemed to us to add more complexity than to clarify the physics. We prefer to stick to the usual definition of the superadiabatic Rayleigh number involving $\Delta T - \Delta T_{sa}$ and, therefore, we must consider the possibility of convection with negative superadiabatic Rayleigh numbers.

### 4.3. Convection simulations

With the same numerical code as in Ricard et al. (2022), we solve the mass, momentum and energy conservation system (16a)-(16b)-(16c) with the software Dedalus (Burns et al., 2020, www.dedalus-project.org), which handles coupled partial differential equations that are

solved using a spectral decomposition for the horizontal direction and Chebyshev polynomials for the vertical direction. The variables are advanced in time by a Runge-Kutta scheme of four-stage and order 3. We choose the number of points so that the top and bottom thermal boundary layers are described by at least 10 Chebyshev collocation nodes. In our simulations, the surface temperature is $T_t = T_0$ and the bottom temperature $T_b = rT_0$. We assume that the surface pressure is $P_t = P_0 = 0$ (hence that the surface density is $\rho_0$).

As discussed in Ricard et al. (2022), it is tricky to work with a compressible fluid on a fixed numerical grid. Indeed, we do not know what the total mass in the convective layer must be to ensure that when convection is well established, the surface pressure is zero. We therefore perform our simulations for given $Ra$, $r$ and $\mathcal{S}$, starting with different initial masses (the initial mass is the integral of an initial adiabatic and hydrostatic profile which itself depends only on the choice of a surface adiabatic density or temperature) until the average surface pressure is statistically zero (the local surface pressure itself remains a function of space and/or time and is commonly interpreted as equivalent to express the presence of a dynamic topography induced by convection).

### 4.4. Heat flux in compressible convection

The efficiency of heat transport by convection is expressed by the Nusselt number, the ratio of the heat flux actually transported to that which would be transported by conduction alone. Since the adiabatic temperature gradient must be established before convection can begin, it is natural to subtract the heat conducted along the adiabat from both the actual heat transport and the conductive heat transport used in the definition of the Nusselt number (e.g., Malkus, 1954; Tilgner, 2011; Jones et al., 2021). A superadiabatic Nusselt number is therefore defined with

$$Nu_{sa} = \frac{Q - Q_a}{\Delta T - \Delta T_a}, \qquad (21)$$

where $Q$ and $Q_a$ are the surface heat flux and the adiabatic surface heat flux (in (21), heat fluxes and temperatures are dimensionless, heat fluxes are normalized by $kT_0/H$ and temperatures by $T_0$). Although this definition has been used for perfect gas EoS, where the adiabatic gradient is uniform, it is also the expression used in solid Earth geophysics when
this gradient is non-uniform (e.g., Bercovici et al., 1989), and the widely used benchmark comparison of compressible anelastic codes by King et al. (2010) also removes the adiabatic gradient in the definition of the Nusselt number. For compressible convection at infinite Pr number, the superadiabatic heat flux and the superadiabatic Rayleigh number are related by

$$Nu_{sa} \propto Ra_{sa}^{1/3}, \qquad (22)$$

in agreement with what is found in the Boussinesq approximation (Malkus, 1954; Grossmann and Lohse, 2001).

The relation (22) is usually proposed in situations where $\Delta T_a$ is smaller and often much smaller than $\Delta T$ which is not necessarily verified when $\mathcal{S}$ is large, as convection can occur even with $\Delta T_a \ge \Delta T$. Notice also that, as shown by the blue dashed line in Fig. 1, the adiabatic heat flow at the surface may be very large and $Q - Q_a$ may be negative even in the case where $\Delta T_a \le \Delta T$. The "adiabatic" Nusselt number, $Nu_a = Q_a/\Delta T_a = \mathcal{T}_a^T/\Delta T_a$ (see (13)) is of order $\mathcal{S}$ as $\mathcal{T}_a^T/\Delta T_a \approx 1$. For $\mathcal{S} = 10$, this is already a large heat flow which requires a Rayleigh number of $10^5 - 10^6$ in the Boussinesq case. Previously we had to consider the case with negative superadiabatic Rayleigh number, now we can have a negative superadiabatic Nusselt number. We could of course alter the definition of the Nusselt number but this would be a cosmetic change that would not change the basic physics.

This situation where the convective heat flow can be less than the adiabatic heat flow, is not unknown and may actually happen in the Earth's core (although inertial, electromagnetic and rotational effects not accounted for in our model become crucial in the core). The conductivity of iron is sufficiently large (Stacey and Loper, 2007; de Koker et al., 2012; Gomi et al., 2013) that in the top part of the core, the heat transported along the adiabat may be larger than that transported by convection (Labrosse et al., 1997; Lister and Buffett, 1998). This would imply the presence of a stratified layer in which the superadiabatic temperature carries the heat flow downward to balance the upward transport along the adiabat. How deep convection interacts with a shallow layer with a large adiabatic gradient can now be discussed using some numerical simulations. A very large number of quantities of interest could be computed from these numerical simulations. Here, we will simply discuss the general pattern of convection and the mean temperature profiles, as well as how the energy is transported throughout the layer.

In the Boussinesq approximation used as a reference model, the heat flow through the fluid is simply

$$Q_{Bo} = \overline{\rho w C_v T} - k \frac{d\overline{T}}{dz}, \qquad (23)$$

where $w$ is the vertical velocity and $T$ the total temperature. The overbar indicates an average of the various quantities horizontally and over time; this heat flow is constant with depth. In the Boussinesq approximation, the density $\rho$ is a constant, just as the heat capacity $C_v$ and the heat capacity at constant volume or constant pressure are not distinguished. In a compressible fluid (see Ricard et al., 2022, for details), the relevant quantity advected by the flow is the enthalpy $\mathcal{H}$ that can be deduced from (4) by integration of $d\mathcal{H} = Td\mathcal{S} + dP/\rho$

$$\mathcal{H} = \left( C_v + \frac{\alpha_0 K_v^0}{\rho} \right) T + \frac{K_v^0}{(n-1)} \frac{\rho^{n-1}}{\rho_0^n}, \qquad (24)$$

and the heat flow can then be written as

$$Q = \overline{\rho w C_v T_{sa}} - k \frac{d\overline{T_{sa}}}{dz} + \overline{\rho w (\mathcal{H} - C_v T_{sa})} - k \frac{dT_a}{dz} - \overline{u \tau_{xz} + w \tau_{zz}}, \qquad (25)$$

where we consider separately the adiabatic, $T_a$, and superadiabatic, $T_{sa} = T - T_a$, temperatures.

The first two terms converge to their counterparts in (23) when the adiabatic temperature is constant, the third term (which can obviously

be simplified with the first) is a correction due to the fact that it is enthalpy rather than specific heat $C_v T$ that is being transported. The fourth term is the conduction along the adiabat and the last is the work flow ($\tau_{xz}$ and $\tau_{zz}$ are the deviatoric stresses, $u$ the horizontal velocity). If we had used dimensionless variables, the ratio of this term to the first would be $\mathcal{S}\Delta T/(\text{Ra}T_0)$ and indeed, would become negligible in the Boussinesq approximation when $\mathcal{S} \to 0$.

### 4.5. Convection at small dissipation number and large temperature drop ($\mathcal{S} = 1, r = 10$)

In a first simulation at a moderate Rayleigh number $Ra_{sa} = 10^8$, we consider a situation in which the adiabatic temperature drop is small compared to the imposed temperature difference, with $\mathcal{S} = 1$ and $r = 10$ which corresponds more or less to terrestrial conditions. A typical snapshot of the superadiabatic temperature field is depicted in Fig. 4. This situation remains relatively close to the classical Rayleigh-Bénard convection in the Boussinesq approximation but some differences are notable. Due to dissipation, the descending and ascending plumes are rather discontinuous. They tend to form clusters, which is consistent with the suggestion by Schubert et al. (2004) that the two superplume regions observed in the deep Earth's mantle may be clusters of smaller plumes which coalesced into a large region of hot, buoyant material.

The time-average temperature profile is depicted in blue in Fig. 5. The dots along the temperature profile correspond to the nodes of the Chebyshev polynomials used by the Dedalus software. Top and bottom boundary layers have comparable thicknesses. The adiabatic profile (red) is quasi linear and provides a close approximation to the real temperature (Fig. 5). The total mass under the actual temperature profile and under the adiabatic profile is the same and leads to a zero average pressure at the surface when the convection is statistically steady.

In this simulation the time averaged heat flux is $Q = 431$, $Q_a = 4.41$, $\Delta T_{sa} = 6.78$ and the Nusselt number, $Nu_{sa} = 62.92$ and therefore $Nu_{sa} = 0.14 Ra_{sa}^{1/3}$. The prefactor is in agreement with other simulations performed in the Boussinesq regime (Sotin and Labrosse, 1999) or in the fully compressible case with an ideal gas Eos (Curbelo et al., 2019). We depict the profiles of the various components of the heat flow in Fig. 6 (see (25)). In the left panel (a), we plot the total heat flow (green), the transport of specific heat (red) and the conduction along the superadiabatic temperature profile (blue) (i.e., the total heat flow $Q$ and the first two terms of eq. (25)). The total heat flow (green) would be depth-independent when averaged over a very long period. The red and blue components (the specific heat transport and the conductive term) are the only terms carrying heat in the Boussinesq approximation (see (23)). Due to compressibility, other minor contributors to the energy transport are present in panel 6b. The transport of $\mathcal{H} - C_v T_{sa}$ (red), the conduction along the adiabat (blue), and the work flow (green) (i.e., the third, fourth and fifth terms of (25)). These different terms tend to increase the energy transport near the surface and decrease it at depth. Compressibility does not affect the energy transport very much as $\mathcal{S}$ remains small compared to the applied total temperature difference across the layer and the minor components of panel b have amplitudes at most of

### Superdiabatic Temperature

<figure><img src="image_6.png" /><figcaption>Fig. 4. Snapshot of the superadiabatic temperature in the convective layer with $Ra_{sa} = 10^8$, $\mathcal{S} = 1$ and $r = 10$. The downwellings and upwellings are less stable and less continuous with depth than in the Boussinesq approximation. The plumes are not homogeneously distributed but tend to form clusters.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

**Fig. 5.** Temperature (blue) and adiabatic (red) profile in a convection model with Ra<sub>sa</sub> = 10<sup>8</sup>, *D* = 1 and *r* = 10. The fluid mass that can be computed from the adiabatic and hydrostatic profile is by construction the total mass inside the convective fluid. The average pressure at the surface and over time is zero. The two boundary layers have similar thicknesses. The temperature overshoot near the hot bottom boundary layer is slightly more pronounced than that under the lithosphere. The total temperature jump across the layer is *r* − 1 = 9 and the adiabatic temperature jump is 2.22. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

≈ 5% of the major components of panel a.

4.6. Convection at large dissipation number and large temperature drop
(*D* = *r* = 10)

We can now consider the case where the dissipation number becomes comparable to the temperature ratio, *D* = *r* = 10. A typical snapshot of the superadiabatic temperature is depicted in Fig. 7. When the dissipation number is large, it is difficult for cold plumes to continuously

traverse the convective layer (Alboussière et al., 2022). As previously observed by Hansen et al. (1993), hot instabilities gain buoyancy as they rise in the mantle as the thermal expansivity increases. They are stronger and more stationary than the cold instabilities that lose buoyancy with depth. The temperature profiles are shown in Fig. 8. Due to the large heat flow carried out along the adiabat, the cold boundary layer (the lithosphere) is poorly defined and its thickness increases compared to the Boussinesq case. On the contrary, the bottom boundary layer is much less affected because the adiabatic gradient is minimal here. The adiabatic profile (red) and the total temperature (blue) have approximately similar curvatures. We recall that both thermal profiles lead to the same total mass in the mantle; the red adiabatic curve has not been computed to provide a best fit to the observed temperature. Furthermore, the idea that the actual temperature profile should be adiabatic is based on the assumption that the dissipation is negligible which is not the case when *D* is large.

In this simulation the heat flow is $Q = 270$, $Q_a = 30.75$, $\Delta T_{sa} = 5.56$ and the Nusselt number, $Nu_{sa} = 43.03$ and therefore $Nu_{sa} = 0.09 Ra_{sa}^{1/3}$. The prefactor of this equation is smaller than usually found: the Nusselt number, which is inversely related to the thickness of the cold boundary layer, is small because this thickness is increased by the large flux carried by conduction along the adiabat.

As in Fig. 6, the various components of the heat flow are depicted in Fig. 9 when the dissipation number is now *D* = 10. The transport of

Superadiabatic Temperature

<figure><img src="image_8.png" /></figure>

**Fig. 7.** Snapshot of the superadiabatic temperature in the convective layer with Ra<sub>sa</sub> = 10<sup>8</sup>, *D* = 10 and *r* = 10. The downwellings are still less stable and less continuous with depth than the hot plumes but for both hot and cold plumes, crossing the mantle becomes difficult.

<figure><img src="image_9.png" /></figure>

**Fig. 6.** Energy flux through the convective layer. In panel (a), we plot the total energy flow, averaged over time, in green. This flux should become independent of depth if we had performed the averaging over an infinitely long time window. The two main contributors of the energy flux are the advection of specific heat (red) and the conduction along the superadiabatic thermal profile (blue). The former is effective in the bulk of the fluid, the latter in the boundary layers. Other minor contributors to energy transport (panel (b)) are conduction along the adiabat (blue), work flow (green), the contribution due to enthalpy rather than specific heat being transported (red). Notice the difference in scale between the two panels. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Temperature profile (blue) and adiabatic (red) in a convection model with Ra = 10⁸, $\mathcal{D}$ = 10 and r = 10. The mass of the fluid that can be computed from the adiabatic profile is by construction the total mass of the fluid. The total average pressure at the surface, average over time is zero. The thickness of the top boundary is strongly affected by the ability of the fluid to carry heat along the adiabat. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

specific heat (red, panel a), as for the previous case, underestimates the energy transport in the upper part of the mantle and overestimates the energy transport in the lower part. The heat conducted along the superadiabatic gradient (blue, panel (a)) is lower across the top cold boundary where significant heat (12% of the surface heat flow) is carried out along the adiabat (blue curve, panel (b)). Energy transfert due to the difference between enthalpy and specific heat (red, panel (b)) and work flow (green, panel (b)) are also significant. These three minor components added to the two major components in panel (a), lead to a global heat flow (green, panel (a)), independent of depth.

### 4.7. Convection at large dissipation number and small temperature drop ($\mathcal{D}$ = 10 and r = 2.5)

As we discussed, convection can also occur with a negligible superadiabatic temperature jump $\Delta T_{sa}$. We therefore perform a simulation, always at Ra$_{sa}$ = 10⁸ but with $\mathcal{D}$ = 10 and r = 2.5. In that case $\Delta T_{sa}$ = 0.30. A temperature snapshot is shown in Fig. 10. Rising plumes are strong and vigorous but few downwelling currents are visible. In fact, no cold boundary exists and on the contrary a significant portion of the surface corresponds to a maximum of superadiabatic temperature. These characteristics are also obvious in Fig. 11 depicting the average temperature profile (blue). The adiabatic temperature (red) has a slightly lower surface temperature than the real surface temperature.

The various components of the energy transport are shown in Fig. 12. The balance between the various terms is very different from the previous cases (see Figs. 6 and 9). In this simulation the heat flow is $Q = 19.1$, $Q_a$ = 10.7, and the superadiabatic Nusselt number as defined in (28) is now very small $Nu_{sa}$ = 28. When the adiabatic gradient has a large curvature with depth, the usual Ra-Nu relation cannot be used. Conduction along the superadiabatic temperature profile (blue, panel b) becomes less important than along the adiabat (blue panel b). The advection of specific heat (red, panel a) transports the energy in the deep convective layer, but the main transport occurs at shallow depth, along the adiabat (blue, panel b). The work flow (green, panel b) and the distinction between enthalpy and specific heat (red, panel b) only

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Snapshot of the superadiabatic temperature in the convective layer ($\mathcal{D}$ = 10 and r = 2.5). The hot rise plumes spread under an even hotter (in terms of superadiabatic temperatures) top boundary layer.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Energy flow across the convective layer. In panel (a), we plot the total energy flow, averaged on time, in green. The other minor contributors of energy transport are the conduction along the adiabat (blue, panel (b)), the work flow (green, panel (b)), the contribution due to the fact the enthalpy transport is not restricted to the $C_vT$ term (red, panel (b)). Notice the difference in scale between the two panels. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Temperature profile (blue) and adiabatic (red) in a convection model with Ra = $10^8$, $\mathcal{D} = 10$ and $r = 2.5$. Except for a bottom boundary layer the superadiabatic average temperature (blue) follows the general curvature of the adiabat (red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

transports a small excess of energy in the upper part of the convective layer.

## 4.8. Convection with negative superadiabatic Rayleigh number ($\mathcal{D} = 5$ and $r = 1.5$)

As the marginal stability analysis suggests, convection can also occur with a negative superadiabatic Rayleigh number. This is the case when $\mathcal{D} = 5$ and $r = 1.5$ which are parameters that belong to the blue shaded area of the phase diagram of Fig. 3. We choose a negative Rayleigh number, $Ra = -10^8$, smaller than the negative critical number computed in Fig. 3. The resulting convection pattern is typically that shown in Fig. 13. The shallow layer now appears as a warm (non adiabatic) stable layer beneath which vanish faint rising plumes, warmer

than the deep average mantle but colder than the stable surface lid. The temperature profile of Fig. 14 (blue) is far away from the adiabat (red). The temperature profile is linear in the stable conductive layer ($0.7 \le z \le 1$) but a weak hot boundary still remains at the bottom. The heat transport components (Fig. 15) are really dominated by the conductive terms along the superadiabatic profile (panel (a)) and along the adiabat (panel (b)). These two terms largely cancel each other out but together extract the heat flow transported in the deep layers by convection at the surface (red, panel (a)). The work flow and the difference between enthalpy and specific heat are totally negligible. In this case the Nusselt number is positive with negative numerator and denominator; $Nu = 12.5$ ($Q = 0.76$, $Q_a = 3.98$, $\Delta T = 0.5$, $\Delta T_a = 0.76$), but with a negative Rayleigh number.

## 5. Conclusions

In this article, we examined how convection is affected by compressibility when the dissipation number, which is a measure of the non-Boussinesq effects, is large. This situation occurs in Super-Earths, i. e. in solid or liquid planets whose radius is much larger than that of the Earth. Such planets are quite common and the surprising large variety of exoplanets that have been found so far suggests that planets even larger than those considered here (say up to 3 times the radius of the Earth) exist.

The characteristics of compressible flows are controlled by an adia-

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Snapshot of the superadiabatic temperature in a convection model with a negative Rayleigh number $Ra = -10^8$, $\mathcal{D} = 5$ and $r = 1.5$. As shown in Fig. 3 convection occurs with these parameters. The colour scale is chosen to emphasize the weak rising plumes and saturates near the surface.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Energy flow across the convective layer. In panel (a), we plot the total energy flow, averaged on time, in green. The conduction along the superadiabatic temperature is totally different from the previous cases as now, it carries heat downward near the surface (blue panel (a)). The advection of specific heat (red, panel (a)) dies out at shallow depth. The other minor contributors of energy transport are the conduction along the adiabat (blue, panel (b)), the work flow (green, panel (b)), the contribution due to the fact the enthalpy transport is not restricted to the $C_pT$ term (red, panel (b)). Notice the difference in scale between the two panels. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Temperature profile (blue) and adiabatic (red) in a convection model with Ra = −10⁸, $\mathcal{D}$ = 5 and r = 1.5. The top to bottom temperature ratio is now significantly smaller than the adiabatic temperature difference. Convection occurs with a negative Rayleigh number. A conductive layer in which the temperature gradient is constant fills the upper 20% of the convective layer but a weak bottom boundary layer allows the emergence of plumes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** The conduction terms along the superadiabatic temperature, in blue panel (a), and along the adiabat, in blue panel (b), largely cancel each other. At depth the energy is transported by convection (red, panel (a)). The other terms of energy transport are negligible. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

batic state that provides an approximate reference to the thermodynamic values during developed convection. In a first section we discussed some properties of the adiabatic conditions when the fluid obeys a simple Murnaghan EoS with a Grüneisen parameter decreasing inversely with density. This EoS is fairly simple but faithfully reproduces the behavior of solids and liquids at high pressure and temperature. Some consequences of this EoS for adiabatic profiles of density or temperature are independent of the nature of the fluid (the properties of solid silicates or liquid iron properties can be reproduced by this generic EoS). This EoS imposes that incompressibility increases rapidly with pressure. The effects of compression are therefore concentrated in the upper layers while the deep layers can hardly be compressed any further. This imposes a strong curvature on adiabatic density and tem-

perature. It is unlikely that the bottom adiabatic temperature is more than twice the adiabatic temperature at the surface, even in large Super-Earths. Although the adiabatic temperature gradient, which plays a major role in convection, can be very large at the surface, it remains moderate at depth and, surprisingly, it decreases, rather than increases, with the dissipation number $\mathcal{D}$. Another way of looking at this is to note that, although the surface dissipation number $\mathcal{D}$ estimated from the radius of a planet can be very large, Murnaghan's EoS requires its mean to be less than the Grüneisen parameter, i.e., typically less than about 1.

We then study the marginal stability of compressible convection. For the EoS we have chosen, the strong curvature of the adiabat facilitates convection in the deep layers and favors conductive heat transport in shallow layers. Convection can develop in deep layers even when the total temperature difference between the bottom and top surfaces is equal to or less than the adiabatic temperature difference. Convection with a negative superadiabatic Rayleigh number is therefore possible.

We then explore different cases of developed convection. Our computations are carried out with a moderate superadiabatic Rayleigh number (|Ra| = 10⁸), various dissipation numbers and various bottom to top temperature ratios. The simulations are carried out for a fluid without inertia which is only valid for the creeping convection of planetary mantles, but without any further approximation. There is, of course, no indication of the appropriate Rayleigh numbers and temperature ratios for Super-Earths. This temperature ratio (and the Rayleigh number) depends on the surface and bottom conditions of the planet's mantle. The former depends on atmospheric composition, distance from the star, and the planet's internal dynamics. The latter depends on the planet's formation mechanisms, core segregation, mantle radioactive content and planet's formation age. It is probably unlikely that the temperature ratio imposed across silicate mantles is as small as in the cases of sections 4.7 or 4.8 ($\mathcal{D} \gg r \approx 1$) but there are maybe more strange situations than are dreamt in our philosophy. The situation discussed in section 4.6 ($\mathcal{D} \approx r \gg 1$) could be a more common situation. In this case, the top and bottom boundary layers are of very different thicknesses, and the concept of a top boundary layer, a lithosphere, may lose its meaning as significant conduction along the adiabat may suppress or dampen convection. The heat transport in the case of large dissipation number can occur through various combinations of enthalpy transport (which cannot be limited to specific heat), conduction along the adiabatic gradient and along the average superadiabatic temperature gradient and a work flow term.

In this article, viscosity has been assumed to be uniform throughout the mantle. The various behaviors that have been described, as a function of temperature ratio *r* and dissipation number $\mathcal{D}$, are solely due to the response of the EoS to thermal forcing. This means that we do not expect the large scale behavior of convection to depend on the exact rheology. Of course, the rheology can be different – non-newtonian and non-uniform – and the precise small-scale flow structures will depend on this, however the overall picture of convection is mainly governed by the EoS. Numerical experiments of compressible convection in Super-Earths with non-uniform viscosity and thermal conductivity profiles have been conducted by Kameyama and Yamamoto (2018). They were carried out using the truncated anelastic liquid approximations that lead to thermodynamic inconsistencies. Solving the exact equations with depth-dependent viscosity and conductivity appears straightforward, however, using Dedalus and very minor modifications to our code.

Although our simulations without inertia are valid only for creeping flows, i.e., for solid state convection, some of our conclusions hold for the case where convection occurs in liquids, including liquid metals. First, as already mentioned, the chosen EoS is a much better starting point for correctly expressing the thermodynamical equations that control the flow than is often done. The characteristics of adiabatic properties discussed in this paper remain valid for fluids. Of course, other terms, inertia, rotation, electromagnetic effects would have to be added. For these cases, anelastic approximations should be used for numerical modeling. It remains necessary to start from a realistic EoS
and carefully check the consistency of the approximations. Situations like in sections 3.7 or 3.8 where the dissipation number becomes much larger than the temperature ratio, are probably common and could prevail in the Earth's core. Note that in these cases, an estimate of the extracted heat flow based on the adiabatic gradient is meaningless, as convection does not develop near the surface (see the difference between the adiabatic transport, blue curve of Fig. 15(b) and the total energy flow, green curve of Fig. 15(a)). The superadiabatic temperature drives the heat flow down (blue curve of Fig. 15(a)) and the surface heat flow remains comparable to what deep convection (red curve of Fig. 15(a)) is capable of transporting.

## Author statement

* Yanick Ricard and Thierry Alboussière conceived the main conceptual ideas.

* Yanick Ricard wrote the numerical code and the original version of this article.

* Thierry Alboussière verified the analytical basis and improved the manuscript by comments and suggestions.

## Declaration of Competing Interest

No conflict of interest.

## Data availability

No data was used for the research described in the article.

## References

Alboussière, T., Ricard, Y., 2013. Reflections on dissipation associated with thermal convection. J. Fluid Mech. 725, 749-751.

Alboussière, T., Ricard, Y., 2017. Rayleigh-benard stability and the validity of quasi-boussinesq or quasi-anelastic liquid approximations. J. Fluid Mech. 817, 264-305.

Alboussière, T., Curbelo, J., Dubuffet, F., Labrosse, S., Ricard, Y., 2022. A playground for compressible natural convection with a nearly uniform density. J. Fluid Mech. 940, A9.

Alfè, D., Price, G., Gillan, M., 2002. Iron under Earth's core conditions: liquid-state thermodynamics and high-pressure melting curve from ab initio calculations. Phys. Rev. B 65 (16).

Anderson, O.L., 1979. Evidence supporting the approximation $γρ$ = const for the Grüneisen parameter of the Earth lower mantle. J. Geophys. Res. 84, 3537-3542.

Bercovici, D., Schubert, G., Glatzmaier, G.A., Zebib, A., 1989. 3-dimensional thermal-convection in a spherical-shell. J. Fluid Mech. 206, 75-104.

Braginsky, S., Roberts, P., 1995. Equations governing convection in earth's core and the geodynamo. Geophys. Astrophys. Fluid Dyn. 79 (1-4), 1-97.

Burns, K.J., Vasil, G.M., Oishi, J.S., Lecoanet, D., Brown, B.P., 2020. Dedalus: a flexible framework for numerical simulations with spectral methods. Phys. Rev. Res. 2, 023068.

Chandrasekhar, S., 1957. An Introduction to the Study of Stellar Structure. Dover Publications, Astrophysical monographs.

Choblet, G., Cadek, O., Couturier, F., Dumoulin, C., 2007. OEDIPUS: a new tool to study the dynamics of planetary interiors. Geophys. J. Int. 170 (1), 9-30.

Curbelo, J., Duarte, L., Alboussière, T., Dubuffet, F., Labrosse, S., Ricard, Y., 2019. Numerical solutions of compressible convection with an infinite Prandtl number: comparison of the anelastic and anelastic liquid models with the exact equations. J. Fluid Mech. 873, 646-687.

de Koker, N., Steinle-Neumann, G., Vlcek, V., 2012. Electrical resistivity and thermal conductivity of liquid Fe alloys at high p and t, and heat flux in earth's core. Proc. Natl. Acad. Sci. U. S. A. 109 (11), 4070-4073.

Gomi, H., Ohta, K., Hirose, K., Labrosse, S., Caracas, R., Verstraete, M.J., Hernlund, J.W., 2013. The high conductivity of iron and thermal evolution of the earth's core. Phys. Earth Planet. Inter. 224, 88-103.

Grossmann, S., Lohse, D., 2001. Thermal convection for large Prandtl numbers. Phys. Rev. Lett. 86, 3316-3319.

Gubbins, D., Alfè, D., Masters, G., Price, G.D., Gillan, M.J., 2003. Can the Earth's dynamo run on heat alone? Geophys. J. Int. 155 (2), 609-622.

Hansen, U., Yuen, D.A., Kroening, S.E., Larsen, T.B., 1993. Dynamic consequences of depth-dependent thermal expansivity and viscosity on mantle circulations and thermal structure. Phys. Earth Planet. Inter. 77, 205-223.

Jarvis, G.T., McKenzie, D.P., 1980. Convection in a compressible fluid with infinite prandtl number. J. Fluid Mech. 96, 515-583.

Jeffreys, H., 1930. The instability of a compressible fluid heated below. Proc. Cambridge Phil. Soc. 26 (2), 170-172.

Jones, C.A., Mizerski, K.A., Kessar, M., 2021. Fully developed anelastic convection with no-slip boundaries. J. Fluid Mech. 930.

Kameyama, M., Yamamoto, M., 2018. Numerical experiments on thermal convection of highly compressible fluids with variable viscosity and thermal conductivity: implications for mantle convection of super-Earths. Phys. Earth Planet. Inter. 274, 23-36.

King, S.D., Lee, C., van Keken, P.E., Leng, W., Zhong, S., Tan, E., Tosi, N., Kameyama, M. C., 2010. A community benchmark for 2-D Cartesian compressible convection in the Earth's mantle. Geophys. J. Int. 180 (1), 73-87.

Labrosse, S., Poirier, J., LeMouel, J., 1997. On cooling of the earth's core. Phys. Earth Planet. Inter. 99 (1-2), 1-17.

Lantz, S., Fan, Y., 1999. Anelastic magnetohydrodynamic equations for modeling solar and stellar convection zones. Astron. J. 121 (1), 247-264.

Leng, W., Zhong, S., 2008. Viscous heating, adiabatic heating and energetic consistency in compressible mantle convection. Geophys. J. Int. 173 (2), 693-702.

Lister, J., Buffett, B., 1998. Stratification of the outer core at the core-mantle boundary. Phys. Earth Planet. Inter. 105 (1-2), 5-19.

Malkus, W., 1954. The heat transport and spectrum of thermal turbulence. Proc. Royal Soc. London Series A-Math. Phys. Sci. 225 (1161), 196-212.

Miyagoshi, T., Kameyama, M., Ogawa, M., 2015. Thermal convection and the convective regime diagram in super-Earths. J. Geophys. Res. Planets 120 (7), 1267-1278.

Murnaghan, F.D., 1951. Finite Deformation of an Elastic Solid. John Wiley and Sons, York.

Murphy, C.A., Jackson, J.M., Sturhahn, W., 2013. Experimental constraints on the thermodynamics and sound velocities of hcp-fe to core pressures. J. Geophys. Res. 118 (5), 1999-2016.

Ogura, Y., Phillips, N., 1962. Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci. 19 (2), 173-179.

Otegi, J.F., Bouchy, F., Helled, R., 2020. Revisited mass-radius relations for exoplanets below 120 m-circle plus. Astron. Astrophys. 634.

Petit, A.T., Dulong, P.L., 1819. Recherches sur quelques points importants de la théorie de la chaleur. Ann. Chim. Phys. 10, 395-413.

Rayleigh, J., 1916. On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. London, Edinburgh, Dublin Philos. Mag. J. Sci. 32 (192), 529-546.

Ricard, Y., 2015. 7.02 - physics of mantle convection. In: Schubert, G. (Ed.), Treatise on Geophysics, second edition. Elsevier, Amsterdam, pp. 31-87.

Ricard, Y., Alboussière, T., Labrosse, S., Curbelo, J., Dubuffet, F., 2022. Fully compressible convection for planetary mantles. Geophys. J. Int. 230 (2), 932-956.

Schubert, G., Turcotte, D., Olson, P., 2001. Mantle Convection in the Earth and Planets. Cambridge University Press.

Schubert, G., Masters, G., Olson, P., Tackley, P., 2004. Superplumes or plume clusters? Phys. Earth Planet. Inter. 146 (1), 147-162. Plumes and Superplumes.

Schwarzschild, K., 1906. Ueber das gleichgewicht der sonnenatmosphärischen Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen. Math. Phys. Klasse 1906, 41-53.

Sotin, C., Labrosse, S., 1999. Three-dimensional thermal convection in an iso-viscous, infinite prandtl number fluid heated from within and from below: applications to the transfer of heat through planetary mantles. Phys. Earth Planet. Inter. 112 (3-4), 171-190.

Stacey, F.D., Davis, P.M., 2004. High pressure equations of state with applications to the lower mantle and core. Phys. Earth Planet. Inter. 142, 137-184.

Stacey, F.D., Loper, D.E., 2007. A revised estimate of the conductivity of iron alloy at high pressure and implications for the core energy balance. Phys. Earth Planet. Inter. 161 (1-2), 13-18.

Tackley, P.J., 2008. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid. Phys. Earth Planet. Inter. 171 (1-4, SI), 7-18.

Tilgner, A., 2011. Convection in an ideal gas at high rayleigh numbers. Phys. Rev. E 84(2, 2).

Zhang, S.X., Yuen, D.A., 1996. Various influences on plumes and dynamics in time-dependent, compressible mantle convection in 3-D spherical shell. Phys. Earth Planet. Inter. 94, 241-267.

Zhong, S., McNamara, A., Tan, E., Moresi, L., Gurnis, M., 2008. A benchmark study on mantle convection in a 3-D spherical shell using CitcomS. Geochem. Geophys. Geosyst. 9.