Onset of thermal convection in a solid spherical shell with melting at either or both boundaries

SUMMARY

Thermal convection in planetary solid (rocky or icy) mantles sometimes occurs adjacent to liquid layers with a phase equilibrium at the boundary. The possibility of a solid–liquid phase change at the boundary has been shown to greatly help convection in the solid layer in spheres and plane layers and a similar study is performed here for a spherical shell with a radius-independent central gravity subject to a destabilizing temperature difference. The solid–liquid phase change is considered as a mechanical boundary condition and applies at either or both horizontal boundaries. The boundary condition is controlled by a phase change number, Φ, that compares the timescale for latent heat exchange in the liquid side to that necessary to build a topography at the boundary. We introduce a numerical tool, available at https://github.com/amorison/stablinrb, to carry out the linear stability analysis of the studied setup as well as other similar situations (Cartesian geometry, arbitrary temperature and viscosity depth-dependent profiles). Decreasing Φ makes the phase change more efficient, which reduces the importance of viscous resistance associated to the boundary and makes the critical Rayleigh number for the onset of convection smaller and the wavelength of the critical mode larger, for all values of the radii ratio, γ. In particular, for a phase change boundary condition at the top or at both boundaries, the mode with a spherical harmonics degree of 1 is always favoured for Φ ≲ 10⁻¹. Such a mode is also favoured for a phase change at the bottom boundary for small (γ ≲ 0.45) or large (γ ≳ 0.75) radii ratio. Such dynamics could help explaining the hemispherical dichotomy observed in the structure of many planetary objects.

Instability analysis, Numerical solutions, Mantle processes

1 INTRODUCTION

In several instances, convection occurs in solid layers of planetary interiors with a phase equilibrium at either or both of the horizontal boundaries. Such a situation allows for the possibility of mass transfer across the boundary which modifies the dynamics and heat transfer characteristics of convection (Labrosse et al. 2018). The first situation in which this question has been addressed concerns convection in the inner core of the Earth (Alboussière et al. 2010; Monnereau et al. 2010; Deguen et al. 2013, 2018; Mizzon & Monnereau 2013). These papers showed that the phase change at the inner-core boundary can be included in the dynamics of the inner core as a phase change boundary condition (BC). This BC combines dynamic pressure, radial velocity and its radial derivative, with a control parameter, the dimensionless phase change number, denoted by Φ in the present paper. This parameter is the ratio between the timescale associated with heat transfer in the liquid and the timescale for building a topography of the boundary by a perturbation of the stress field in the solid. For a large value of Φ, the topography of the boundary is limited by its buoyancy and the radial velocity effectively decays to zero, which is the classical no-penetration boundary condition. On the other end, for a small value of Φ, the boundary topography is suppressed by melting and freezing associated with latent heat transfer on the liquid side, which allows matter to flow through by changing phase. These papers showed that the critical value for the onset of convection decreases linearly with Φ at low values of this parameter and that the preferred mode of convection in that case has a spherical harmonic degree of one, without deformation, which corresponds to a translation of the inner core.

In other planetary situations, a solid spherical shell can be the seat of convection with the possibility of melting and freezing at either or both of its boundaries. The ice layers of satellites of Jupiter and Saturn are generally bounded by a liquid water ocean, below for the surface ice layers or above for the high pressure (HP) ice layers in the large satellites, Titan and Ganymede (Hussmann et al. 2015). The implications of the phase change boundary condition on convection in HP ice layers have been studied by use of numerical models by Lebec et al. (2023a) and Lebec et al. (2023b). These papers show that the efficiency of heat and mass transfer by convection in these ice layers is increased by the presence of a solid–liquid phase equilibrium at the top. Moreover, they show that the thermal structure is quite different from that with a classical no-penetration BC: only one boundary layer exists, on the no-penetration boundary, so that only hot plumes exist in that case, with a broad down-welling return flow. They also show that, in some cases, melting of the HP ice layer could occur at the bottom, in contact with the solid core, as was also shown by previous models that did not include the phase change BC (Choblet et al. 2017; Kalousová et al. 2018; Kalousová & Sotin 2018). In the case of a very salty ocean, a dense ocean could stabilize between the HP ice layer and the rocky mantle (Lebec et al. 2023b) which would require to apply a phase change BC at both boundaries.

During their early stages, rocky planets like Earth are thought to undergo large amounts of melting, in particular associated with giant impacts, which leads to their magma ocean stage (see Elkins-Tanton 2012, for a review). In most cases, the crystallization is thought to proceed upward from the bottom and the question of the onset of convection in the solid requires to consider the effect of the phase change boundary condition (Morison et al. 2019). The persistence of a dense molten layer at the bottom of Earth mantle termed the basal magma ocean, has been suggested to explain the present state as imaged by seismology (Labrosse et al. 2007). Such a layer has recently been proposed to exist today at the bottom of the Martian mantle (Samuel et al. 2023; Khan et al. 2023). In that case, a phase change BC needs also to be considered at the bottom when modelling convection in the solid mantle. Moreover, extrapolating the present situation to the early times of such planets leads to the possibility of a solid mantle bounded above and below by magma oceans (Labrosse et al. 2007). Irrespective of the mechanism by which such a situation can emerge during formation of the planet (Labrosse et al. 2015), the dynamics of the solid layer needs to include the phase change BC at both boundaries in that case (Bolrão et al. 2021).

Some effort has been done to extend the work done in the context of inner core dynamics on the effect of the phase change BC on convection. The simplest situation is undoubtedly the plane layer case. Labrosse et al. (2018) studied that case using linear and weakly non-linear analysis and showed phase change helps to start convection and leads to large scale flow. When both boundaries are subject to a phase change, a translation mode of convection is possible and the deforming mode has a very long wavelength. These predictions have been tested and extended by Agrusta et al. (2019) using a mantle convection model.

The spherical shell geometry has been considered by Deguen (2013) who studied the linear onset of convection in a solid spherical shell with solid–liquid phase change at either or both of its horizontal boundaries. Deguen (2013) confirmed the important effect of such boundary condition on the onset: compared to the classical non-penetrative BC, allowing mass transfer across the boundary makes the critical Rayleigh number smaller and the preferred mode at onset larger scale (smaller spherical harmonic degree). Deguen (2013) considered a gravity acceleration that varies linearly with radial direction, as applies to the Earth core. This allowed carrying out a linear stability analysis in a fully analytical fashion. However, in the present study, we consider a planetary mantle with constant gravity. The linear stability analysis can then be solved via numerical tools as presented in this paper.

In this paper, we only study the solid part of the crystallizing mantle. The presence of magma oceans is parametrized with boundary conditions at the top and bottom boundaries of the solid domain. Such a vision is of course overly simplistic, but it allows us to isolate and understand interesting consequences of the presence of (magma) oceans on solid-state convection. Several previous papers have considered the effect of variable compositions associated with fractional crystallization at the boundaries (Morison et al. 2019; Bolrão et al. 2021; Lebec et al. 2023b) but we focus here on thermal convection only, using linear stability analysis.

2 FORMULATION OF THE PROBLEM

The solid part of the crystallizing mantle is considered to be a spherical shell of external radius R⁺ and internal radius R⁻. Its thickness is denoted L = R⁺ − R⁻. The geometry is considered constant, leaving aside the possibility of net freezing or melting that has been treated by Morison et al. (2019). For generality purpose, we consider global magma oceans can exist above and/or below the solid layer. They are also considered to be spherical shells. The Top Ocean (TO) extends from the external solid boundary R⁺ to the surface of the planet at radius R_T; the interface between the TO and the solid layer is called the Top Ocean–Mantle Boundary (TOMB). The Basal Ocean (BO) extends from the Core–Mantle Boundary (CMB) at radius R_c to the internal solid boundary R⁻; the interface between the BO and the solid layer is called the Bottom Ocean–Mantle Boundary (BOMB). See Fig. 1 for a visual representation of these layers.

Figure 1.

Equatorial view of the problem setup. We consider a solid layer (green) surrounded by global magma oceans (red) above and/or below that solid. The radii associated with the various interface positions are labelling the relevant arrows.

This section presents the conservation equations describing mantle convection and the boundary conditions to apply to them, with a focus on the effect of the solid–liquid phase change. Finally, making these equations and associated boundary conditions dimensionless exhibits the controlling parameters of this simplified system.

2.1 Conservation equations

The velocity u = (u_r, u_θ, u_ϕ), pressure P and temperature T fields in the solid are linked together by a set of conservation equations. For the sake of simplicity, the viscosity η and thermal diffusivity κ are assumed to be uniform throughout the mantle. The gravity acceleration |$\mathbf {g}=-g\hat{\mathbf {r}}$| is radial and its intensity g is uniform. This work is done under the Boussinesq approximation where all density variations are considered negligible except, of course, in the buoyancy term which is the motor of thermal convection. Moreover, the solid mantle has an extremely large Prandtl number, which means inertia terms are negligible compared to the viscous forces. We hence consider the solid mantle is an infinite Prandtl number fluid with no inertia. Under these assumptions, the mass, momentum and energy conservation equations are the following (e.g. Chandrasekhar 1961):

$$\begin{eqnarray} \nabla \cdot \mathbf {u} = 0, \end{eqnarray}$$

(1)

$$\begin{eqnarray} \mathbf {0} = -\boldsymbol{\nabla }P + \eta \nabla ^2\mathbf {u} - \rho _\alpha g\hat{\mathbf {r}}, \end{eqnarray}$$

(2)

$$\begin{eqnarray} \dfrac{\partial ^{}{T}}{\partial {t}^{}} + \mathbf {u}\cdot \boldsymbol{\nabla }T = \kappa \nabla ^2 T, \end{eqnarray}$$

(3)

where ρ_α is the density used in the buoyancy term. It is related to the temperature field via the thermal expansion coefficient α which is considered constant. Denoting ρ the reference density at a temperature T₀, the buoyancy density ρ_α can be expanded as

$$\begin{eqnarray} \rho _\alpha = \rho \left(1 - \alpha (T-T_0)\right). \end{eqnarray}$$

(4)

Note moreover that in this problem, we do not consider internal heating occurring in the solid. Only heat advection and diffusion are included in eq. (3). This is justified here by the fact that heat producing elements (U, K, Th) are incompatible (e.g. Huang et al. 2013) and therefore are expected to reside in the magma oceans rather than the solid.

2.2 Boundary conditions

The boundary positions r = R⁻ (BOMB) and r = R⁺ (TOMB) correspond to the case where no convection operates in the solid, and the solid and liquid layers are at thermodynamic equilibrium. In this case—referred to as ‘the motionless position’ in this paper—both boundaries are spherical, and at the phase change temperatures T⁻ and T⁺. Moreover, no heat is gained or lost by the overall system (which would result in a net melting or freezing of the magma oceans). In practice, if convection operates in the solid layer, matter at the boundary departs from its motionless position and therefore forms a topography h with respect to it. This topography is either a solid residing in a liquid or a liquid residing in a solid of the same composition, and therefore prone to melting or freezing. This phase change can act as an erosion mechanism of the topography, effectively resulting in an exchange of matter between the solid and the liquid layers. Such a system comprises a large wealth of complexities and therefore several assumptions need be made in order to ease its study. Fig. 2 illustrates what happens as the actual solid–liquid boundary is deformed around its motionless position. The topography with respect to the motionless position is denoted h⁺ at the TOMB and h⁻ at the BOMB; it is by convention positive when oriented towards a higher radial position. The actual liquid/solid interface of the system is then at r = R⁻ + h⁻ for the bottom boundary and r = R⁺ + h⁺ for the top boundary.

Figure 2.

Dynamic topography at the boundaries. Note that the curvature has been removed and the height of topography exaggerated for readability purpose. The horizontal dot-dashed lines at z = ±1/2 correspond to the motionless position of the boundaries r = R^± when no convection operates in the solid and there is no net freezing or melting of the magma oceans. Matter departing from the boundary creates a topography h. The thick line on the right side is the super-isentropic temperature profile in the motionless steady state. The temperature at the topography follows the melting temperature (thin solid line), and therefore departs from the isentropic profile in the liquid. Figure from Labrosse et al. (2018).

The normal stress continuity at either boundary can be written as:

$$\begin{eqnarray} 2\eta \dfrac{\partial ^{}{u_r}}{\partial {r}^{}}(h_s) - P(h_s) = -P(h_l), \end{eqnarray}$$

(5)

where f(h_s) is the quantity f at the solid side of the boundary, and f(h_l) is the quantity f at the liquid side of the boundary. Note that the viscosity of the liquid is several orders of magnitude smaller than that of the solid; viscous forces on the liquid side are therefore neglected. The pressure is then written as the sum between the hydrostatic pressure |$\bar{P}$| and the dynamic pressure p:

$$\begin{eqnarray} &&2\eta \dfrac{\partial ^{}{u_r}}{\partial {r}^{}}(h_s) - p(h_s) \overbrace{-\bar{P}_0+\rho _s gh}^{-\bar{P}(h_s)}\\ &&\quad=- p(h_l) \underbrace{-\bar{P}_0+\rho _l gh}_{-\bar{P}(h_l)}, \end{eqnarray}$$

(6)

where |$\bar{P}_0$| is the hydrostatic pressure at the boundary in the motionless position. Moreover, the dynamic pressure p(h_l) in the liquid can be neglected since its time average value is null and the dynamics of the liquid part is expected to be fast compared to the one of the solid. Eq. (6) then becomes

$$\begin{eqnarray} 2\eta \dfrac{\partial ^{}{u_r}}{\partial {r}^{}}(h_s) - p(h_s) +\Delta \rho gh = 0 \end{eqnarray}$$

(7)

where Δρ = ρ_s − ρ_l is the density contrast between the solid and the liquid. Note that for the system to be mechanically stable, Δρ should be positive at the top boundary and negative at the bottom boundary. The last variable that needs to be connected to the dynamics of the solid mantle is the topography h. At the TOMB, a positive (resp. negative) topography h⁺ is formed by solid (liquid) matter that goes toward the liquid (solid) part but that has not melted (crystallized) yet. Conversely, at the BOMB, a positive topography h⁻ corresponds to liquid matter going toward the solid but that has not frozen yet.

The topography h^± is formed by the velocity of the solid itself u_r as well as the freezing-front velocity V_r. This freezing-front velocity is chosen positive along the outward radial direction. A positive value of V_r corresponds to freezing at the TOMB, and melting at the BOMB. The Lagrangian derivative of h^± hence verifies

$$\begin{eqnarray} \dfrac{\mathrm{D}h^\pm }{\mathrm{D}t} = u_r + V_r. \end{eqnarray}$$

(8)

Moreover, energy conservation at the boundary is given by Stefan’s law:

$$\begin{eqnarray} \rho _sL_hV_r = (\mathbf {q}_l - \mathbf {q}_s)\cdot \hat{\mathbf {r}} \end{eqnarray}$$

(9)

where L_h is the latent heat of crystallization and q is the heat flux on either side of the boundary. Assuming the topography is at thermodynamic equilibrium (crystallization on a pre-existing freezing front with a moderate freezing speed argues against the importance of undercooling), the temperature at the boundary T(h) is the melting temperature at the relevant pressure. It is related to the melting temperature at the motionless position of the boundary T^± as follows:

$$\begin{eqnarray} T(h^\pm ) = T^\pm - \dfrac{\partial ^{}{T_L}}{\partial {P}^{}}\rho _lgh^\pm . \end{eqnarray}$$

(10)

|$\dfrac{\partial ^{}{T_L}}{\partial {P}^{}}$| is the Clapeyron slope (T_L being the phase change temperature), considered constant on the pressure range across the topography h.

Solving this set of equations at both boundaries requires a full convection model in the liquid to determine the heat flux in the liquid q_l. However, the timescale at which the liquid evolves is much shorter than that of the solid. This makes solving Navier–Stokes in a consistent way for both the solid and liquid layers impractical with our current computing power. The liquid behaviour at the timescale at which solid-state convection operates is therefore parametrized. The parametrization we use is the one that was introduced for the core by Deguen et al. (2013). Assuming advection dominates heat transport, one can write

$$\begin{eqnarray} \mathbf {q}_l\cdot \hat{\mathbf {r}}\sim \rho _l C_{p,l}v^{\prime }\delta T \end{eqnarray}$$

(11)

where v′ is the typical fluid velocity in the liquid and δT is the temperature departure from the average temperature profile in the liquid at the topography. The average temperature profile is assumed to be isentropic as convection in magma oceans is expected to be vigorous (e.g. Solomatov 2015). This leads to

$$\begin{eqnarray} \delta T = \left( \left.\dfrac{\partial ^{}{T}}{\partial {P}^{}}\right|_S - \dfrac{\partial ^{}{T_L}}{\partial {P}^{}}\right)\rho _lgh. \end{eqnarray}$$

(12)

Moreover, the heat flux from the solid |$\mathbf {q}_s\cdot \hat{\mathbf {r}}$| is considered small compared to the one carried out in the liquid and is therefore neglected in eq. (9). Injecting eqs (11) and (12) in eq. (9) gives the following expression for the freezing velocity V_r:

$$\begin{eqnarray} V_r = - \dfrac{h}{\tau _\phi } \end{eqnarray}$$

(13)

with τ_ϕ the phase change timescale

$$\begin{eqnarray} \tau _\phi =\dfrac{\rho _s L_h}{\rho _l^2gC_{P,l}v^{\prime }(\partial _PT_L-\left.\partial _PT\right|_S)}. \end{eqnarray}$$

(14)

τ_ϕ is the timescale at which heat is carried in the ocean from places of freezing to places of melting.

Eq. (8) can be simplified by realizing what are the various timescales at play. Let U be the typical scale for the velocity in the solid. Denoting τ_η = η/(ΔρgL) the viscous timescale, eq. (7) gives us Uτ_η as scale for the topography. Let us introduce the timescale τ_c associated with the convective derivative of the topography in eq. (8), which is the minimum between two values: the timescale associated to changes in convection in the solid and l_h/U, the ratio of the horizontal length scale and the typical solid flow velocity. Using these scales to render eq. (8) dimensionless yields:

$$\begin{eqnarray} \frac{\tau _\eta }{\tau _c} \frac{\mathrm{D}h^{\prime \pm }}{\mathrm{D}t^\prime } = u_r^\prime - \frac{\tau _\eta }{\tau _\phi } h^\prime , \end{eqnarray}$$

(15)

primed variables being dimensionless. Following the example of Earth’s mantle for which the timescale of mantle convection is much larger than that for postglacial rebound, we can safely assume that τ_c ≫ τ_η, meaning the time derivative of topography in eq. (8) can be neglected (Dh/Dt ∼ 0). Plugging eq. (13) in eq. (8) yields the following expression for the topography h

$$\begin{eqnarray} h = u_r\tau _\phi . \end{eqnarray}$$

(16)

Substituting h with its expression in eq. (7) and assuming the topography is small compared to the thickness of the domain, one obtains the following boundary condition applied at r = R^±

$$\begin{eqnarray} \Delta \rho g \tau _\phi u_r + 2\eta \dfrac{\partial ^{}{u_r}}{\partial {r}^{}} - p = 0. \end{eqnarray}$$

(17)

The topography h is now an implicit variable of the problem that can be computed a posteriori. Note that such a boundary condition allows a degree-1 topography to develop (the so-called translation modes discussed further in this paper), and therefore can lead to a displacement of the centre of mass. As detailed by Deguen et al. (2013), gravitational potential perturbations can be incorporated in the pressure leading to the same formulation when considering self-gravity.

Assuming viscous forces in the liquid are negligible, the shear stress continuity at each boundary is written:

$$\begin{eqnarray} \tau _{r\theta } &=& \eta \left(r\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({\dfrac{u_\theta }{r}}\right) +\dfrac{1}{r}\dfrac{\partial ^{}{u_r}}{\partial {\theta }^{}} \right) = 0,\\ \tau _{r\phi } &=& \eta \left(r\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({\dfrac{u_\phi }{r}}\right) +\dfrac{1}{r\sin \theta }\dfrac{\partial ^{}{u_r}}{\partial {\phi }^{}} \right) = 0. \end{eqnarray}$$

(18)

Those equations are applied directly to the reference boundary R^±, the topography h being considered small.

Finally, boundary conditions are required for the temperature field. As shown by Labrosse et al. (2018), in the limit of a small topography, we can assume that the solid–liquid phase transition occurs at a laterally constant temperature. Hence, Dirichlet conditions are used at both boundaries:

$$\begin{eqnarray} T(r=R^-) = T^-, \quad T(r=R^+) = T^+. \end{eqnarray}$$

(19)

2.3 Dynamic pressure choice

Applying the boundary condition eq. (17) at both boundaries might lead to a non-zero average radial velocity at a boundary. This would be equivalent to a net-freezing of one magma ocean and a net-melting of the other. We want to decouple completely net-freezing or melting of the magma oceans (related to the long term evolution of the system) and dynamic freezing or melting at the boundaries (related to the dynamic topography formed by viscous forces). To this effect, the dynamic pressure needs to be chosen to ensure the boundary condition eq. (17) cannot yield a non-zero average radial velocity. In the rest of this document, 〈•〉 denotes the lateral mean of a quantity. As previously, + and − superscripts denote quantities evaluated at the top and bottom boundary, respectively.

Denoting T_p the temperature profile such that T = T_p yields p = 0, the Stokes eq. (2) becomes

$$\begin{eqnarray} \mathbf {0} = -\boldsymbol{\nabla }p + \eta \nabla ^2\mathbf {u} + \rho g\alpha (T-T_p)\hat{\mathbf {r}}. \end{eqnarray}$$

(20)

Integrating the Stokes equation over the entire solid domain Ω and projecting on |$\hat{\mathbf {r}}$| gives:

$$\begin{eqnarray} 0 &=& \int _+ \left( 2\eta \dfrac{\partial ^{}{u_r}}{\partial {r}^{}} - p \right){\mathrm{d}}S - \int _- \left( 2\eta \dfrac{\partial ^{}{u_r}}{\partial {r}^{}} - p \right){\mathrm{d}}S \\ &&+\,\, \int _\Omega \rho g\alpha (T-T_p){\mathrm{d}}V. \end{eqnarray}$$

(21)

Plugging the phase change boundary condition eq. (17) in eq. (21) leads to

$$\begin{eqnarray} &&\Delta \rho ^+g\tau _\phi ^+(R^+)^2\left\langle {u_r} \right\rangle ^+ + |\Delta \rho ^-|g\tau _\phi ^-(R^-)^2\left\langle {u_r} \right\rangle ^-\\ &&\quad= \rho g\alpha \int _{R^-}^{R^+}\left(\left\langle {T} \right\rangle -T_p\right)r^2{\mathrm{d}}r. \end{eqnarray}$$

(22)

Finally, mass conservation gives us:

$$\begin{eqnarray} (R^+)^2\left\langle {u_r} \right\rangle ^+ = (R^-)^2\left\langle {u_r} \right\rangle ^-. \end{eqnarray}$$

(23)

This leads to

$$\begin{eqnarray} \left( \Delta \rho ^+\tau _\phi ^+ + |\Delta \rho ^-|\tau _\phi ^- \right)(R^+)^2\left\langle {u_r} \right\rangle ^+ = \rho \alpha \int _{R^-}^{R^+}\left(\left\langle {T} \right\rangle -T_p\right)r^2{\mathrm{d}}r. \\ \end{eqnarray}$$

(24)

Eq. (24) shows the average topography at the boundaries |$\tau _\phi ^\pm \left\langle {u_r} \right\rangle ^\pm$| is directly proportional to the average buoyancy of the bulk. Choosing T_p equal to 〈T〉 therefore ensures the average topography (and average radial velocity) is zero at all times. Defining the dynamic pressure as p = P − 〈P〉, the Stokes eq. (2) becomes

$$\begin{eqnarray} \mathbf {0}=-\boldsymbol{\nabla }p + \eta \nabla ^2\mathbf {u} + \rho g\alpha (T-\left\langle {T} \right\rangle )\hat{\mathbf {r}}. \end{eqnarray}$$

(25)

2.4 Dimensionless equations

The equations are made dimensionless in order to reduce the number of parameters describing the physical problem. The scales for distance, time, and temperature are respectively the thickness of the domain L = R⁺ − R⁻, the thermal diffusion timescale |$\frac{L^2}{\kappa }$|⁠, and the temperature difference ΔT = T⁻ − T⁺ between the two interfaces. The dimensionless temperature |$\widetilde{T}$| is defined as:

$$\begin{eqnarray} \widetilde{T} = \dfrac{T - T^+}{\Delta T} \end{eqnarray}$$

(26)

so that it is always between 0 and 1.

Using the same symbols for dimensionless and dimensional quantities, the non-dimensional conservation equations are

$$\begin{eqnarray} \nabla \cdot \mathbf {u}=\mathbf {0}, \end{eqnarray}$$

(27)

$$\begin{eqnarray} \mathbf {0} = -\boldsymbol{\nabla }p + \nabla ^2\mathbf {u} + \mathrm{Ra}(T-\left\langle {T} \right\rangle )\mathbf {\hat{r}}, \end{eqnarray}$$

(28)

$$\begin{eqnarray} \dfrac{\partial ^{}{T}}{\partial {t}^{}} + \mathbf {u}\cdot \boldsymbol{\nabla }T = \nabla ^2 T, \end{eqnarray}$$

(29)

with Ra the Rayleigh number defined as:

$$\begin{eqnarray} \mathrm{Ra}\equiv \dfrac{\rho g\alpha \Delta TL^3}{\eta \kappa }. \end{eqnarray}$$

(30)

This dimensionless number compares the buoyancy forces which drive convection to the momentum and heat diffusion coefficients which hinder it.

With eq. (26), the boundary conditions for the temperature are straightforward:

$$\begin{eqnarray} T^- &=& T(R^-) = 1; \\ T^+ &=& T(R^+) = 0. \end{eqnarray}$$

(31)

The free-slip boundary condition eq. (18) gives:

$$\begin{eqnarray} r\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({\dfrac{u_\theta }{r}}\right) +\dfrac{1}{r}\dfrac{\partial ^{}{u_r}}{\partial {\theta }^{}} &=& 0,\\ r\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({\dfrac{u_\phi }{r}}\right) +\dfrac{1}{r\sin \theta }\dfrac{\partial ^{}{u_r}}{\partial {\phi }^{}} &=& 0. \end{eqnarray}$$

(32)

Finally, the normal stress continuity condition eq. (17) leads to:

$$\begin{eqnarray} \pm \Phi ^{\pm }u_r + 2\dfrac{\partial ^{}{u_r}}{\partial {r}^{}} - p = 0. \end{eqnarray}$$

(33)

Φ⁺ and Φ⁻ are the dimensionless phase change numbers of the top and bottom interfaces, defined as

$$\begin{eqnarray} \Phi ^{\pm } \equiv \dfrac{|\Delta \rho |^\pm g L}{\eta }\tau _\phi . \end{eqnarray}$$

(34)

These parameters are ratios between the phase change timescale τ_ϕ and the timescale needed to build topography by viscous forces. These numbers represent the resistance of the melting/freezing boundaries to flow of matter through them. If Φ → 0 at one interface, the boundary condition eq. (33) reduces to a normal stress free condition, meaning the interface is fully permeable. Physically, melting and freezing of matter is much quicker than viscous building of topography, allowing matter to pass easily through the interface. The height of the topography is indeed limited by the rate at which it is melted/frozen away instead of its buoyancy. On the contrary, if Φ → ∞, the boundary condition eq. (33) imposes that u_r = 0 and becomes equivalent to the definition of dynamic topography under the usual non-penetrative boundary condition applied in convection models (see Ricard et al. 2014, for a discussion). Physically, viscous building of topography is fast enough for the height of the topography to be limited by its weight rather than melting or freezing, hence preventing matter from crossing the boundary.

After making the equations dimensionless, only four parameters are necessary to describe the system. The inner radius R⁻ characterizing the geometry (the outer radius being R⁺ = 1 + R⁻ since lengths are made dimensionless with the thickness of the domain), the Rayleigh number Ra describing the strength of convection, and the two phase change numbers Φ^± parametrizing the behaviour of the two interfaces with the magma oceans.

3 METHOD

A first approach to study the system described previously is a linear stability analysis. We study here the stability of the conductive state in the solid, varying the phase change numbers Φ^± as well as the geometry of the system defined by the dimensionless position of the internal boundary R⁻. The approach presented here is akin to that used by Deguen (2013) who studies a similar problem with three differences:

1. The gravity acceleration is considered constant in our study while it varies linearly with the radial position in Deguen (2013);

2. We neglect internal heating in our study, while Deguen (2013) considers a volumetrically heated domain;

3. The Rayleigh number is defined with the temperature difference and the domain thickness in our study, but with the volumetric heating rate and the outer radius in Deguen (2013); comparison of the critical Rayleigh number values between the two studies should therefore be made with care.

The choices made by Deguen (2013) were dictated by the will to use the same approach as that used for the inner core for which they fully apply. The choices done here are more relevant to convection in solid mantles and require a different solution technique, as shown below. The numerical tool presented in this paper is nonetheless very versatile and can treat the cases solved by Deguen (2013) with minor modifications.

3.1 Motionless reference state

The system of partial differential equations formed by the conservation eqs (27)–(29) and their boundary conditions eqs (32) and  (31) admits a purely conductive (i.e. motionless) steady solution. This solution, denoted by an overline, is defined by

$$\begin{eqnarray} \bar{\mathbf {u}} = \mathbf {0} \end{eqnarray}$$

(35)

and

$$\begin{eqnarray} \nabla ^2\bar{T} = 0. \end{eqnarray}$$

(36)

The latter equation along with the boundary conditions on temperature eq. (31) takes as solution

$$\begin{eqnarray} \bar{T} = \dfrac{R^+R^-}{r}-R^-. \end{eqnarray}$$

(37)

Injecting this in the Stokes eq. (28) gives

$$\begin{eqnarray} \boldsymbol{\nabla }\bar{p} = \mathbf {0}. \end{eqnarray}$$

(38)

Therefore, |$\bar{p}$| is constant throughout the solid domain. Since it should be null to satisfy the phase change boundary condition eq. (33), one obtains

$$\begin{eqnarray} \bar{p} = 0 . \end{eqnarray}$$

(39)

Introducing the temperature anomaly |$\Theta = T-\bar{T}$|⁠, the conservation equations become:

$$\begin{eqnarray} \nabla \cdot \mathbf {u} = 0, \end{eqnarray}$$

(40)

$$\begin{eqnarray} \mathbf {0} = -\boldsymbol{\nabla }p + \nabla ^2\mathbf {u} + \mathrm{Ra}\left(\Theta -\left\langle {\Theta } \right\rangle \right)\mathbf {\hat{r}}, \end{eqnarray}$$

(41)

$$\begin{eqnarray} \dfrac{\partial ^{}{\Theta }}{\partial {t}^{}} + \mathbf {u}\cdot \boldsymbol{\nabla }(\Theta +\bar{T}) = \nabla ^2\Theta . \end{eqnarray}$$

(42)

The boundary conditions on Θ are Θ = 0 at both boundaries.

3.2 Poloidal potential formulation

Since the fluid is considered isoviscous and incompressible, the velocity field can be reduced to a scalar field, the poloidal potential |$\mathcal {P}$| defined as (e.g. Ricard & Vigny 1989; Ribe 2015):

$$\begin{eqnarray} \mathbf {u} = \boldsymbol{\nabla }\times \boldsymbol{\nabla }\times (\mathbf {r}\mathcal {P}). \end{eqnarray}$$

(43)

One can notice that the poloidal potential is related to the stream vector |$\boldsymbol{\Psi }$|⁠:

$$\begin{eqnarray} \mathbf {u}=\boldsymbol{\nabla }\times \boldsymbol{\Psi }\iff \boldsymbol{\Psi }=\boldsymbol{\nabla }\times (\mathbf {r}\mathcal {P}). \end{eqnarray}$$

(44)

The three components of the velocity field are then:

$$\begin{eqnarray} u_r &=& \dfrac{1}{r} \mathcal {L}^2\mathcal {P}, \end{eqnarray}$$

(45)

$$\begin{eqnarray} u_\theta &=& \dfrac{\partial ^{}{}}{\partial {\theta }^{}}\left({\dfrac{1}{r}\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({r\mathcal {P}}\right)}\right), \end{eqnarray}$$

(46)

$$\begin{eqnarray} u_\phi &=& \dfrac{1}{\sin \theta }\dfrac{\partial ^{}{}}{\partial {\phi }^{}}\left({\dfrac{1}{r}\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({r\mathcal {P}}\right)}\right). \end{eqnarray}$$

(47)

|$\mathcal {L}^2$| is the scalar operator defined as:

$$\begin{eqnarray} \mathcal {L}^2\bullet &\equiv& \dfrac{\partial ^{}{}}{\partial {r}^{}}\left({r^2\dfrac{\partial ^{}{\bullet }}{\partial {r}^{}}}\right) - r^2\nabla ^2\bullet \\ &=& -\dfrac{1}{\sin \theta }\dfrac{\partial ^{}{}}{\partial {\theta }^{}}\left({\sin \theta \dfrac{\partial ^{}{\bullet }}{\partial {\theta }^{}}}\right) -\dfrac{1}{\sin ^2\theta }\dfrac{\partial ^{2}{\bullet }}{\partial {\phi }^{2}}. \end{eqnarray}$$

(48)

Using the following properties (e.g. Dormy 1997):

$$\begin{eqnarray} \mathcal {L}^2\bullet = (\boldsymbol{\nabla }\times \boldsymbol{\nabla }\times (\bullet \mathbf {r}))\cdot \mathbf {r}, \end{eqnarray}$$

(49)

$$\begin{eqnarray} \boldsymbol{\nabla }\times \boldsymbol{\nabla }\times \mathbf {u} = -\boldsymbol{\nabla }\times \boldsymbol{\nabla }\times (\mathbf {r}\nabla ^2\mathcal {P}), \end{eqnarray}$$

(50)

one can find the poloidal formulation of the momentum conservation equation by taking twice the curl of eq. (41):

$$\begin{eqnarray} \mathcal {L}^2\left(\nabla ^4\mathcal {P}- \dfrac{\mathrm{Ra}}{r}\Theta \right) = 0. \end{eqnarray}$$

(51)

Assuming the poloidal potential and temperature fields can be developed with spherical harmonics (see Section 3.3), solutions to eq. (51) involving a non-zero |$(\nabla ^4\mathcal {P}- \frac{\mathrm{Ra}}{r}\Theta )$| are purely radial functions of harmonic degree 0. By definition, Θ does not have a degree-0 component. Degree-0 components in |$\mathcal {P}$| can be taken as zero without loss of generality since they have no role in eq. (43), leaving

$$\begin{eqnarray} \nabla ^4\mathcal {P}= \dfrac{\mathrm{Ra}}{r}\Theta . \end{eqnarray}$$

(52)

Introducing |$\mathcal {Q}$| such as:

$$\begin{eqnarray} \mathcal {Q}\equiv \nabla ^2\mathcal {P}, \end{eqnarray}$$

(53)

eq. (52) can be written

$$\begin{eqnarray} \nabla ^2\mathcal {Q}= \dfrac{\mathrm{Ra}}{r}\Theta . \end{eqnarray}$$

(54)

The heat eq. (42) becomes:

$$\begin{eqnarray} \dfrac{\partial ^{}{\Theta }}{\partial {t}^{}}+(\boldsymbol{\nabla }\times \boldsymbol{\nabla }\times (\mathbf {r}\mathcal {P}))\cdot \boldsymbol{\nabla }\Theta = \nabla ^2\Theta + R^-R^+\dfrac{1}{r^3}\mathcal {L}^2\mathcal {P}. \end{eqnarray}$$

(55)

Using eqs (45) to (47), the free-slip boundary condition eq. (32) leads to:

$$\begin{eqnarray} \dfrac{\partial ^{2}{\mathcal {P}}}{\partial {r}^{2}} + (\mathcal {L}^2- 2)\dfrac{\mathcal {P}}{r^2} = C^\pm \end{eqnarray}$$

(56)

where C^± denotes an arbitrary constant. The choice of C^± does not matter to perform the linear stability analysis since it is a term of harmonic degree 0 and therefore vanishes when equations are written for higher harmonic degrees. Degree-0 terms are ignored since they are forbidden by our definition of the dynamic pressure as shown in Section 2.3.

Finally, since the dynamic pressure does not appear in the poloidal formulation of the momentum conservation eq. (54), it should be eliminated from the normal stress continuity boundary condition eq. (33). Denoting |$\boldsymbol{\omega }=\boldsymbol{\nabla }\times \mathbf {u}$| the vorticity and recalling ∇ · u = 0, one can deduce that

$$\begin{eqnarray} \nabla ^2\mathbf {u} = -\boldsymbol{\nabla }\times \boldsymbol{\omega }. \end{eqnarray}$$

(57)

The projection of the momentum conservation eq. (41) along |$\boldsymbol{\hat{\theta }}$| gives

$$\begin{eqnarray} \dfrac{1}{r}\dfrac{\partial ^{}{p}}{\partial {\theta }^{}} + (\boldsymbol{\nabla }\times \boldsymbol{\omega })\cdot \mathbf {\boldsymbol{\hat{\theta }}} = 0. \end{eqnarray}$$

(58)

Since there is no source of toroidal potential in the studied problem, there is no radial vorticity (e.g. Ribe 2015). Hence,

$$\begin{eqnarray} (\boldsymbol{\nabla }\times \boldsymbol{\omega })\cdot \mathbf {\boldsymbol{\hat{\theta }}} = -\dfrac{1}{r}\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({r\omega _\phi }\right). \end{eqnarray}$$

(59)

With eqs (45) and (46), one obtains:

$$\begin{eqnarray} \omega _\phi = (\boldsymbol{\nabla }\times \mathbf {u})\cdot {\boldsymbol {\hat{\phi }}}= \dfrac{1}{r}\left(\dfrac{\partial ^{}{}}{\partial {r}^{}}\left({ru_\theta }\right) - \dfrac{\partial ^{}{u_r}}{\partial {\theta }^{}}\right)= \dfrac{\partial ^{}{}}{\partial {\theta }^{}}\left({\nabla ^2\mathcal {P}}\right). \end{eqnarray}$$

(60)

Putting eqs (59) and (60) in eq. (58) leads to:

$$\begin{eqnarray} \dfrac{\partial ^{}{p}}{\partial {\theta }^{}} = \dfrac{\partial ^2}{\partial \theta \partial r} \left(r\nabla ^2\mathcal {P}\right). \end{eqnarray}$$

(61)

Performing a similar calculation on the ϕ direction gives this relation between the dynamic pressure and the poloidal potential:

$$\begin{eqnarray} -p + \dfrac{\partial ^{}{}}{\partial {r}^{}}\left({r\nabla ^2\mathcal {P}}\right) = K(r) \end{eqnarray}$$

(62)

where K(r) is a purely radial function constant along the θ and ϕ directions. This allows us to substitute the pressure in the phase change condition eq. (33):

$$\begin{eqnarray} \pm \Phi ^{\pm }\dfrac{1}{r}\mathcal {L}^2\mathcal {P}+ \dfrac{\partial ^{}{}}{\partial {r}^{}}\left({\dfrac{2}{r}\mathcal {L}^2\mathcal {P}-r\mathcal {Q}}\right) = K^\pm . \end{eqnarray}$$

(63)

Note that as C^± in eq. (56), the actual value of K^± is irrelevant in this study.

3.3 Spherical harmonics development

Perturbations of the poloidal potential |$\mathcal {P}$|⁠, |$\mathcal {Q}$| and the temperature field Θ are developed using spherical harmonics as following:

$$\begin{eqnarray} \mathcal {P}&=& \sum _{l=1}^{\infty }\sum _{m=-l}^l P_l^m(r) Y_l^m(\theta ,\phi )e^{\sigma _l t}, \\ \mathcal {Q}&=& \sum _{l=1}^{\infty }\sum _{m=-l}^l Q_l^m(r) Y_l^m(\theta ,\phi )e^{\sigma _l t}, \\ \Theta &=& \sum _{l=1}^{\infty }\sum _{m=-l}^l T_l^m(r) Y_l^m(\theta ,\phi )e^{\sigma _l t}. \end{eqnarray}$$

(64)

Note that the l = 0 harmonic is not taken into account since it corresponds to the motionless conductive state.

In the frame of linear stability analysis, one can study each mode (l, m) in an independent way. For a given problem (R⁻, Φ⁺, Φ⁻), the goal of the analysis is to determine which mode is the most unstable and what is its associated critical Rayleigh number Ra_c. Moreover, the problem is degenerated in terms of lateral orientation. Hence, the growth rate σ_l of any given mode (l, m) only depends on l. For readability purposes, m indices are dropped.

The differential operators can easily be written for any mode l. Indeed, spherical harmonics are eigenfunctions of the |$\mathcal {L}^2$| operator. Applying |$\mathcal {L}^2$| to a given mode reduces to:

$$\begin{eqnarray} \mathcal {L}^2\bullet = l(l+1)\bullet . \end{eqnarray}$$

(65)

The Laplacian operator applied to a given mode is then:

$$\begin{eqnarray} \mathcal {D}^{2}_l\bullet \equiv \dfrac{\partial ^{2}{\bullet }}{\partial {r}^{2}}+\dfrac{2}{r}\dfrac{\partial ^{}{\bullet }}{\partial {r}^{}}- \dfrac{l(l+1)}{r^2}\bullet . \end{eqnarray}$$

(66)

The conservation equations eqs (53) to (55) can be written as (the only neglected non-linear term being the advection |$\mathbf {u}\cdot \boldsymbol{\nabla }\Theta$|⁠):

$$\begin{eqnarray} Q_l = \mathcal {D}^{2}_lP_l, \end{eqnarray}$$

(67)

$$\begin{eqnarray} \mathcal {D}^{2}_lQ_l = \dfrac{\mathrm{Ra}}{r}T_l, \end{eqnarray}$$

(68)

$$\begin{eqnarray} (\sigma _l-\mathcal {D}^{2}_l)T_l = R^-R^+\dfrac{l(l+1)}{r^3}P_l. \end{eqnarray}$$

(69)

Finally, the boundary conditions eqs (31), (56) and (63) can be written as:

$$\begin{eqnarray} (T_l)^\pm =0, \end{eqnarray}$$

(70)

$$\begin{eqnarray} \left(\dfrac{{\mathrm{d}}^{2}{P_l}}{{\mathrm{d}}{r}^{2}}\right)^\pm + [l(l+1)-2]\left(\dfrac{P_l}{r^2}\right)^\pm = 0, \end{eqnarray}$$

(71)

$$\begin{eqnarray} \pm \Phi ^\pm l(l+1)\left(\dfrac{P_l}{r}\right)^\pm + \dfrac{{\mathrm{d}}^{}{}}{{\mathrm{d}}{r}^{}}\left({\dfrac{2l(l+1)}{r}P_l-rQ_l}\right)^\pm = 0. \end{eqnarray}$$

(72)

Note that the lateral constants C^± and K^± do not appear in those equations since they are of degree l = 0. Also, as mentioned previously, the linearized equations are independent of the order m of the considered mode (in other words, the pole chosen to define the spherical harmonics has no physical meaning for the problem at hand).

3.4 Eigenvalue formulation

Using a Chebyshev-collocation approach (e.g. Canuto et al. 1985; Guo et al. 2012; Labrosse et al. 2018), the system defined by eqs (67)–(72) can be formulated as a generalized eigenvalue problem. Chebyshev polynomials are used to expand the perturbations along the radial direction. Each vertical mode P_l, Q_l and T_l is entirely characterized by N + 1 Chebyshev–Gauss–Lobatto nodal points at |$z_i = \cos \frac{i\pi }{N}$| with i = 0…N. To map the z ∈ [ − 1, 1] space to the r ∈ [R⁻, R⁺] space, we use the following transformation:

$$\begin{eqnarray} r_i = \dfrac{z_i+1}{2}+R^-. \end{eqnarray}$$

(73)

Each vertical mode can then be represented by a vector with N + 1 components. For example, the toroidal potential vertical mode P_l is represented by P, the vector such as P_i = P_l(r_i). Similarly, Q_l and T_l are represented by Q and T.

With such a formalism, the successive radial derivatives of each vertical mode at any nodal point can be computed with the help of a differentiation matrix d:

$$\begin{eqnarray} \dfrac{\partial ^{k}{\mathbf {P}}}{\partial {r}^{k}} = \mathbf {d}^k \mathbf {P}. \end{eqnarray}$$

(74)

Note that for numerical precision reasons, the powers of d are computed separately instead of directly as the successive powers of d. The differentiation matrices d^k are calculated with the help of a Python adaptation of DMSUITE (Weideman & Reddy 2000). The Python package is available at https://github.com/labrosse/dmsuite.

Denoting |$\underline{\mathbf {r}}$| the diagonal matrix

$$\begin{eqnarray} \underline{\mathbf {r}}\equiv {\begin{pmatrix}r_0 & \quad & \quad \mathbf {0}\\ & \quad \ddots & \quad \\ \mathbf {0} & \quad & \quad r_N \end{pmatrix}}, \end{eqnarray}$$

(75)

the operator |$\mathcal {D}^{2}_l$| can be written as the matrix D²:

$$\begin{eqnarray} \mathbf {D}^2 \equiv \mathbf {d}^2 + 2\underline{\mathbf {r}}^{-1}\mathbf {d} - l(l+1)\underline{\mathbf {r}}^{-2}. \end{eqnarray}$$

(76)

The system defined by eqs (67)–(72) is then equivalent to the matrix equation

$$\begin{eqnarray} \mathbf {L} \mathbf {X} = \sigma \mathbf {R} \mathbf {X} \end{eqnarray}$$

(77)

with

where 1 is the identity matrix. The extra row and column on top and right of the matrices are respectively the column and row indices of each of the submatrices. For example, the top left submatrix of the matrix L is only the first row (hence the 0 on the extra column) of the matrix |$\mathbf {d}^2 + [l(l+1)-2]\underline{\mathbf {r}}^{-2}$|⁠. Similarly, the bottom right submatrix of the matrix L is the matrix D² without its first and last rows and columns. Note that the boundaries of the temperature vertical mode are excluded because the Dirichlet boundary condition eq. (70) is then naturally enforced.

Determining the modes X satisfying eq. (77) as well as the associated growth rate σ is a generalized eigenvalue problem. Given a physical problem Π defined by the three parameters Π = (R⁻, Φ⁺, Φ⁻), and any Rayleigh number Ra and harmonic degree l, the finite eigenvalue σ with the greatest real part is the growth rate of the perturbation of degree l. The eigenvector associated with σ is the vertical modes of the perturbation. For a given physical problem Π and an harmonic l, one can compute the growth rate of the perturbation as a function of the Rayleigh number. The neutral Rayleigh number Ra_n is the Rayleigh number such as |$\mathcal {R}(\sigma (\mathrm{Ra}_n))=0$| where |$\mathcal {R}$| denotes the real part. Finally, for a given problem Π, one can compute the neutral Rayleigh number as a function of the harmonic degree l of the perturbation. The degree l_c for which Ra_n(l) is minimal is the most unstable mode of the problem Π. The associated Rayleigh number Ra_c = Ra_n(l_c) is the critical Rayleigh number for Π.

A numerical implementation of this approach has been developed for this study and is freely available at https://github.com/amorison/stablinrb. This tool is more general than the case presented here as it also allows solving convection in a plane layer (Labrosse et al. 2018), it includes the possibility of net freezing at the boundary (Morison et al. 2019) and permits to consider a vertically varying viscosity and any temperature reference profile, in addition to the conductive one presented above.

4 RESULTS

A simple way to test the linear stability analysis is to perform the analysis with classical no-penetration free-slip condition at both boundaries for which linear analysis results have already been obtained. For comparison, in Cartesian geometry, one gets a critical Rayleigh number Ra_c = 27π⁴/4 ∼ 657.51 and an associated wavenumber |$k_x=\pi / \sqrt{2}$| (Rayleigh 1916). In a spherical shell of aspect ratio γ ≡ R⁻/R⁺ = 0.55, the critical Rayleigh number is Ra_c = 711.95, associated with the harmonic degree l = 3. This is in perfect agreement with the existing literature (e.g. Bercovici et al. 1988). Moreover, when γ → 1, the geometry of the spherical shell tends towards a laterally infinite Cartesian space. With γ = 0.99, the obtained critical Rayleigh number is Ra_c = 657.528 which is very close to the Cartesian value, as expected.

In all the performed analyses, the growth rate σ has no imaginary part. All the modes at Ra_c presented here are hence non-oscillating solutions. This has in fact been mathematically demonstrated for plane layers (Chandrasekhar 1961; Labrosse et al. 2018) but not for the specific case discussed here.

Figs 3 and 4 show the most unstable mode and associated critical Rayleigh number for various cases. Since the linear problem we solve here is degenerated in terms of lateral orientation, all modes with the same harmonic degree l but different orders m have the same growth rate. We choose in Figs 3 and 4 to represent the critical mode of order m = l_c (l_c being the degree of the critical mode) in the equatorial plane |$\theta =\dfrac{\pi }{2}$|⁠. Moreover, instead of showing the poloidal potential perturbation, we show instead isovalues of the stream vector component along |$\boldsymbol{\hat{\theta }}$|⁠. These are streamlines of the flow in the equatorial plane. Enforcing l = m leads to the following relationship between the Ψ_θ component of the stream vector and the poloidal potential eigenmode P_l:

$$\begin{eqnarray} \Psi _\theta (r, \phi ) = \mathcal {R}\left(ilP_l(r)Y_l^l\left( \dfrac{\pi }{2}, \phi \right)\right). \end{eqnarray}$$

(81)

Figure 3.

Temperature perturbation and streamlines of the most unstable mode for no-penetration boundary conditions (left-hand panel) and phase change boundary at the bottom only (right-hand panel). Each row is a different value of the aspect ratio γ of the shell. The critical Rayleigh number Ra_c and harmonic degree l_c are indicated under each figure.

Figure 4.

Temperature perturbation and streamlines of the most unstable mode for a phase change boundary at the top only (left-hand panel) and at both boundaries (right-hand panel). Each row is a different value of the aspect ratio γ of the shell. The critical Rayleigh number Ra_c and harmonic degree l_c are indicated under each figure.

The left column of Fig. 3 shows the classical case with no-penetration boundaries. The critical mode consists in nearly aspect-ratio-one convective rolls, and therefore a critical harmonic degree increasing as the aspect ratio of the shell increases. We can see on the right column of Fig. 3 the effects of a small phase change number at the bottom boundary while a classic non-penetrative condition is prescribed at the top boundary. Owing to the flow-through boundary condition, streamlines are not deviated by the boundary and can instead pass through it as matter melts (downwellings) or freezes (upwellings). Hence, as discussed in Section 2.4, a low value of the phase change number Φ allows matter to cross the boundary. The return current necessary to conserve mass happens in the underlying ocean. This leads to convective patterns with a larger wavelength than the classical rolls, as well as a lower critical Rayleigh number as less deformation is involved in these modes than in the classical rolls.

The left-hand column of Fig. 4 shows that with a flow-through boundary condition with a small phase change number at the top boundary, the most unstable mode is a degree-one pattern regardless of the aspect ratio of the shell. Matter crystallizes on one hemisphere of the shell, goes through the other side avoiding the core, and finally melts on the other hemisphere. As shown in the right-hand column of Fig. 4, when both phase change numbers are small (i.e. both boundaries are flow through), streamlines go straight through the entire shell. All the cases shown on Fig. 4 correspond to degree-one translation modes of convection. The solid shell departs from its motionless position but is constantly recycled and kept in place as it freezes on one side and melts on the other. When both boundaries are flow-through, this translation operates without any deformation in the solid; convection is then only limited by the rate at which melting and freezing can occur. In this last scenario, the critical Rayleigh number is proportional to Φ^±, and can even be arbitrarily small as Φ^± decrease (see Section 5). Finally, when only the top boundary is flow-through, the translation is associated to some deformation in the solid, necessary for the convecting matter to go around the core.

Figs 5 to 7 show the effects of varying the phase change number on the critical Rayleigh number and associated critical harmonic degree. A first observation from these figures is that when both phase change numbers are high (Φ^± ≳ 10²), the system exhibits the same critical Rayleigh numbers and wavelengths as with the classical non-penetrative free-slip boundary conditions. This is expected from the definition of the phase change number Φ as discussed in Section 2.4: for large values of the phase change number, the boundary condition on normal stress eq. (33) converges towards the non-penetrative case. As either or both phase change numbers decrease, the corresponding boundaries transition to the flow-through regime. The critical Rayleigh number and associated harmonic degree decrease as the boundary condition allows for larger wavelength modes of convection. When only one boundary is flow-through (Figs 5 and 6), the critical Rayleigh number decreases of roughly one order of magnitude compared to the no-penetration case with the same geometry. As shown on Fig. 7 and discussed in Section 5, the critical Rayleigh number can be arbitrarily small when both boundaries are flow-through and the translation regime is the most unstable (see Fig. 7). A remarkable feature visible on Figs 6 and 7 is that a small value of Φ⁺ (smaller than about 10) leads to a degree-one translation mode regardless of the aspect ratio of the shell γ or the value of Φ⁻ (a large value for the latter leads to a translation mode with deformation as discussed above).

Figure 5.

Critical Rayleigh number and associated harmonic degree for varying Φ⁻ and various aspect ratios. Left-hand panel: top boundary is non-penetrative, right-hand panel: top boundary is flow-through with Φ⁺ = 10⁻².

Figure 6.

Critical Rayleigh number and associated harmonic degree for varying Φ⁺ and various aspect ratios. Left-hand panel: bottom boundary is non-penetrative, right-hand panel: bottom boundary is flow-through with Φ⁻ = 10⁻².

Figure 7.

Critical Rayleigh number and associated harmonic degree for varying Φ⁺ = Φ⁻ and various aspect ratios.

Finally, Figs 8–10 show the effect of varying the aspect ratio of the shell on the stability of several harmonics degree in three setups: flow-through condition only at the bottom (Fig. 8), only at the top (Fig. 9), and at both boundaries (Fig. 10). As observed before, Figs 9 and 10 show the degree-one mode is the most unstable for any aspect ratio. Moreover, one can notice the neutral Rayleigh number of other modes is much higher than that of the degree-one, showing the degree-one translation mode is strongly favoured. The case of a flow-through boundary only at the bottom depicted in Fig. 8 exhibits an interesting behaviour. For aspect ratio lower than about 0.75, the most unstable mode corresponds to convective rolls that cross the boundary. These rolls are about twice as wide as the classical convective rolls obtained with non-penetrative boundary conditions, as observed in Fig. 3 for γ < 0.8 (for γ = 0.2 it is geometrically impossible for the rolls in the flow-through case to be wider than those in the non-penetrative case since the latter is already of degree one). This is similar to the situation obtained for plane layers (Labrosse et al. 2018) except for the quantization stemming from the sphericity. Increasing γ therefore leads to an increase of l_c. However, for higher aspect ratios of the shell, instead of higher harmonic degrees corresponding merely to these wider rolls, the most unstable mode is of degree-one. This mode is the one shown for γ = 0.8 on Fig. 3, it corresponds to a mode where matter freezes on one inner hemisphere and melts on the other, akin to what happens in the translation mode excepts more deformation is involved in the solid. This new mode has no equivalent in the plane layer situation since it corresponds to an ever increasing wavelength when increasing γ. It takes over from higher-degree modes for γ ≳ 0.75 because it requires less deformation than higher order modes since the flow consists mostly of a rotation around the inner sphere, which is possible because we consider free-slip boundary conditions at both boundaries. Note that this degree-one mode has a critical Rayleigh number close to that of other modes as visible on Fig. 8, allowing competition between the degree-one mode and higher-degree modes at intermediate values of the aspect ratio γ. This contrasts strongly with the translation mode exhibited by cases with a flow-through condition at the top boundary (Figs 9 and 10) that is clearly the most unstable mode. A similar systematics was observed by Deguen (2013) in the case with a gravity acceleration varying linearly with radial position (see his figs 3 and 5).

Figure 8.

Neutral Rayleigh number of several modes as a function of the aspect ratio of the spherical shell. Boundary conditions are flow-through at the bottom (Φ⁻ = 10⁻²) and non-penetrative at the top. Each colour represents a different harmonic degree. At a given aspect aspect ratio, the harmonic with the lowest neutral Rayleigh number is shown as solid, and other harmonics are dotted. In other words, the solid envelope shows the critical Rayleigh number of the system.

Figure 9.

Neutral Rayleigh number of several modes as a function of the aspect ratio of the spherical shell. Boundary conditions are non-penetrative at the bottom and flow-through at the top (Φ⁺ = 10⁻²). Each colour represents a different harmonic degree. The l = 1 harmonic has the lowest neutral Rayleigh number and is shown with a solid line, while other harmonics are shown with dotted lines.

Figure 10.

Neutral Rayleigh number of several modes as a function of the aspect ratio of the spherical shell. Boundary conditions are flow-through at both boundaries (Φ⁺ = Φ⁻ = 10⁻²). Each colour represents a different harmonic degree. The l = 1 harmonic has the lowest neutral Rayleigh number and is shown with a solid line, while other harmonics are shown with dotted lines.

5 ANALYTICAL STUDY OF THE TRANSLATION MODE

The degree-one translation solution is the dominating convection regime at onset when both phase change numbers Φ^± are small. This section presents an analytical determination of the critical Rayleigh number of this mode of convection.

In this section, |$U\hat{\mathbf {t}}$| denotes the translation velocity. U is its amplitude and |$\hat{\mathbf {t}}$| a unit vector along the translation direction. The colatitude is chosen as |$\theta =(\hat{\mathbf {t}},\hat{\mathbf {r}})$|⁠, see Fig. 11 for a schematic of the setup.

Figure 11.

Physical setup and chosen frame in a degree-one translation case. This figure is drawn in the φ = 0 plane. The system is axisymmetric around |$\mathbf {\hat{t}}$| (i.e. φ-invariant).

5.1 Equilibrium between buoyancy and topographic weight

Integrating the Stokes eq. (28) over the entire domain Ω yields

$$\begin{eqnarray} \mathbf {0} &=& \int _+ \left( -p\hat{\mathbf {r}} + \left(\boldsymbol{\nabla }u + \boldsymbol{\nabla }u^T\right)\cdot \hat{\mathbf {r}}\right){\mathrm{d}}S\\ &&-\,\,\int _- \left( -p\hat{\mathbf {r}} + \left(\boldsymbol{\nabla }u + \boldsymbol{\nabla }u^T\right)\cdot \hat{\mathbf {r}}\right){\mathrm{d}}S\\ &&+\,\,\int _\Omega \mathrm{Ra}\left( T-\left\langle {T} \right\rangle \right)\hat{\mathbf {r}}{\mathrm{d}}V. \end{eqnarray}$$

(82)

Plugging the free-slip boundary condition eq. (32) and the phase change boundary condition eq. (33) in eq. (82) leads to

$$\begin{eqnarray} \int _\Omega \mathrm{Ra}\left( T-\left\langle {T} \right\rangle \right)\hat{\mathbf {r}}{\mathrm{d}}V = \int _+\Phi ^+u_r\hat{\mathbf {r}}{\mathrm{d}}S + \int _-\Phi ^-u_r\hat{\mathbf {r}}{\mathrm{d}}S. \end{eqnarray}$$

(83)

Eq. (83) shows the buoyancy of the domain (left-hand side) compensates the total weight of the dynamic topographies (right-hand side).

One can expand the temperature and velocity fields as series of Legendre polynomials. Projecting eq. (83) along the translation direction |$\hat{\mathbf {t}}$| gives the following relation between their degree one components:

$$\begin{eqnarray} &&\mathrm{Ra}\int _{R^-}^{R^+} \underbrace{\int _0^\pi T\cos \theta \sin \theta \mathrm{d}\theta }_{\frac{2}{3}T_1} r^2\mathrm{d}r\\ &&\quad= \left((R^+)^2\Phi ^++(R^-)^2\Phi ^-\right) \underbrace{\int _0^\pi u_r\cos \theta \sin \theta \mathrm{d}\theta } _{\frac{2}{3}{u_r}_1=\frac{2}{3}U}. \end{eqnarray}$$

(84)

This yields the following relation between the translation velocity U and the degree-one component T₁ of the temperature field

$$\begin{eqnarray} U=\dfrac{\mathrm{Ra}}{(R^+)^2\Phi ^++(R^-)^2\Phi ^-} \int _{R^-}^{R^+}T_1r^2\mathrm{d}r. \end{eqnarray}$$

(85)

5.2 Critical Rayleigh number

From eq. (42), the linear growth rate of T₁ is

$$\begin{eqnarray} \dfrac{\partial ^{}{T_1}}{\partial {t}^{}} = -U\dfrac{\partial ^{}{\bar{T}}}{\partial {r}^{}} + \mathcal {D}^2_1T_1 \end{eqnarray}$$

(86)

where |$\bar{T}$| is the reference conductive profile eq. (37), which is also the degree-zero component of the temperature field. The advection of the degree-two component T₂ is neglected as it is a non-linear term. If the Rayleigh number is at the critical value for the degree-one mode of convection (i.e. the translation mode), then the growth rate of T₁ should be zero, since we have seen that it is real. Solving eq. (86) for |$\dfrac{\partial ^{}{T_1}}{\partial {t}^{}}=0$| (with |$T_1^\pm =0$|⁠) and injecting the solution in eq. (85) therefore gives us an equation for the critical Rayleigh number of the translation mode.

$$\begin{eqnarray} &&\mathcal {D}_1^2 T_1 = U\dfrac{\partial ^{}{\bar{T}}}{\partial {r}^{}} = -UR^-R^+\dfrac{1}{r^2}\\ &&\quad\iff T_1 = \dfrac{U}{2}R^-R^+\\ &&\qquad\times \left( 1-\dfrac{(R^+)^2-(R^-)^2}{(R^+)^3-(R^-)^3}r+\dfrac{(R^+)^2(R^-)^3-(R^+)^3(R^-)^2}{(R^+)^3-(R^-)^3}\dfrac{1}{r^2} \right) \end{eqnarray}$$

(87)

This leads to:

$$\begin{eqnarray} \mathrm{Ra}_c = 24\dfrac{\left((R^+)^2\Phi ^++(R^-)^2\Phi ^-\right)\left((R^+)^3-(R^-)^3\right)}{R^+R^-\left((R^+)^2+4R^+R^-+(R^-)^2\right)}. \end{eqnarray}$$

(88)

Introducing the aspect ratio |$\gamma \equiv \dfrac{R^-}{R^+}$|⁠, one obtains:

$$\begin{eqnarray} \mathrm{Ra}_c = 24\dfrac{(\Phi ^++\gamma ^2\Phi ^-)(1-\gamma ^3)}{\gamma (1-\gamma )(\gamma ^2+4\gamma +1)}. \end{eqnarray}$$

(89)

Note that the critical Rayleigh number of the translation mode is directly proportional to Φ⁺ + γ²Φ⁻, and hence to the surface area-weighted average of the phase change numbers. It can therefore be arbitrarily small as the values of these phase-change parameters decrease. Indeed, since no deformation occurs in the solid, the only limiting factor for convection to happen in the translation regime is the rate at which melting and freezing can occur. Moreover, as γ gets close to 1, the geometry of the system gets closer to that of an infinite Cartesian layer. One can notice that

$$\begin{eqnarray} \lim _{\gamma \rightarrow 1}\mathrm{Ra}_c=12(\Phi ^++\Phi ^-), \end{eqnarray}$$

(90)

which is the linear critical Rayleigh number of the translation mode in Cartesian geometry (Labrosse et al. 2018).

6 CONCLUDING REMARKS

Even though a simple approach, the linear stability analysis presented in this paper sheds light on the dramatic consequences of the phase change boundary condition on convection in a solid layer when it is bounded by liquid above, below or both. Note that the results obtained for our system are similar to those of Deguen (2013) for a shell with a gravity increasing with radius and internal heating. This is not surprising given the strong effects of the phase change boundary condition on the behaviour of the flow. These boundary conditions drastically decrease the critical Rayleigh number, affecting the onset of convection in a primitive mantle crystallizing from global magma oceans, as extensively discussed in Morison et al. (2019). The geometry of the flow in the solid is also greatly affected by the possibility for matter to cross the boundary: the expected convective patterns exhibit a much larger wavelength than with classical non-penetrative boundary conditions. In particular, the l = 1 mode is always the most unstable when the top boundary is flow-through, and also for γ ≲ 0.45 or γ ≳ 0.75 when only the bottom boundary is flow-through. In some fashion, this echoes the findings of Zhong & Zuber (2001) in a Mars-like setup: they found that the most unstable mode for Rayleigh–Taylor instability is l = 1 when a rheologically weak asthenosphere is present. This weak layer could play a similar role to a flow-through boundary from the point of view of the rest of the mantle.

Obviously, these results are only strictly valid for the onset of convection since we neglect here the non-linear terms. Zhong et al. (2000); Zhong & Zuber (2001) are other examples of linear stability analyses successfully predicting the most unstable modes of fully non-linear simulations, albeit for different physical setups than the one explored here. Moreover, previous studies using direct numerical simulations of the whole equations set (Agrusta et al. 2019; Bolrão et al. 2021; Lebec et al. 2023a, b) show that the main effects observed in the linear analysis help to understand the fully non-linear results. Indeed, as predicted by the linear stability analysis, small values of the phase change number lead to flow-through boundaries associated with longer wavelength than the typical non-penetrative boundary condition. In Cartesian geometry with phase change at both boundaries, non-linear simulations (Agrusta et al. 2019) present a behaviour that can very well be predicted with a weakly non-linear approach (Labrosse et al. 2018). In spherical geometry, numerical simulations at high Rayleigh number exhibit the translation mode when both boundaries are flow-through (Bolrão et al. 2021).

An aspect that deserves investigation via non-linear simulation is the transition between a flow-through boundary (low Φ) and non-penetrative boundary (high Φ); while the linear stability analysis performed in this study shows that the transition occur over a range of values 10⁻¹ ≲ Φ ≲ 10³, simulations in Cartesian geometry seem to show that the transition in the non-linear regime depends on the ratio Ra/Φ (Agrusta et al. 2019). A systematic exploration of the non-linear dynamics has already been presented in Lebec et al. (2023a) but with a thermal boundary condition at the bottom different from the one considered here, and with only one phase change boundary. Extending these results to the setup of the present paper will be the topic of future work.

ACKNOWLEDGEMENTS

Discussions with Shin-ichi Takehiro were greatly appreciated during the preparation of this paper. The authors thank an anonymous reviewer and Shijie Zhong for their valuable comments which helped improve the clarity of this paper. This study was funded by the Agence Nationale de la Recherche under the grant number ANR-15-CE31-0018-01, MaCoMaOc. AM acknowledges the support of ERC through grant no. 787361-COBOM.

DATA AVAILABILITY

The linear stability code developed for and used in this study is freely available at https://github.com/amorison/stablinrb. It uses the differentiation matrix library freely available at https://github.com/labrosse/dmsuite.