1 Introduction
Rayleigh–Bénard convection is one of the main heat transfer mechanisms in natural sciences, responsible for most of the dynamics of the atmosphere and oceans (Pedlosky Reference Pedlosky1987), plate tectonics (Schubert, Turcotte & Olson Reference Schubert, Turcotte and Olson2001), and dynamo action in planetary cores (Roberts & King Reference Roberts and King2013). It is also one of the most generic examples of pattern formation mechanism in fluid dynamics (e.g. Cross & Hohenberg Reference Cross and Hohenberg1993; Manneville Reference Manneville2004) and has therefore attracted much attention for a century since the work of Lord Rayleigh (Rayleigh Reference Rayleigh1916). However, the mathematical and experimental studies of Rayleigh–Bénard convection have usually considered boundary conditions that are not fully relevant to the natural systems that justified them, their horizontal surfaces being generally considered as subjected to no-slip or free-slip boundary conditions. The former is valid for convection experiments in a tank and for natural fluids bounded by much more viscous envelopes, such as the liquid cores of terrestrial planets and the bottom of the ocean. The latter is often considered as an approximation for a free-surface condition, as applies to a fluid bounded by a much less viscous one. This is in particular the case of the solid planetary mantles that, on long time scales, behave like very viscous fluids (e.g. Turcotte & Oxburgh Reference Turcotte and Oxburgh1967; McKenzie, Roberts & Weiss Reference McKenzie, Roberts and Weiss1974; Jarvis & McKenzie Reference Jarvis and McKenzie1980) and are bounded below and above by liquid or gaseous layers. This approximation neglects the effect of the topography on convection, and some studies have been devoted to the modelling of these effects, which can be dramatic when it is associated with, for example, intense volcanism in hot planets (Monnereau & Dubuffet Reference Monnereau and Dubuffet2002; Ricard, Labrosse & Dubuffet Reference Ricard, Labrosse and Dubuffet2014).
In the present paper, we consider the effects of having horizontal boundaries at which a solid–liquid phase change occurs on Rayleigh–Bénard convection in the creeping solid, that has an infinite Prandtl number (Schubert et al.
Reference Schubert, Turcotte and Olson2001). For simplicity, we consider a Newtonian fluid with a uniform high viscosity, neglecting the effects of more complex rheologies (e.g. Parmentier Reference Parmentier1978; Christensen & Yuen Reference Christensen and Yuen1989; Davaille & Jaupart Reference Davaille and Jaupart1993; Tackley Reference Tackley2000; Bercovici & Ricard Reference Bercovici and Ricard2014), that is bounded by a low-viscosity liquid of the same composition as the convecting solid. The boundary between the liquid and the solid consists of a phase change whose position is controlled by a Clapeyron diagram relating pressure and temperature for phase equilibrium. In the context of planetary interiors, the pressure is largely dominated by the hydrostatic contribution and the interface is on average a horizontal surface. The stress field and associated dynamic pressure due to the dynamics of the solid leads to deformation of the interface with a viscous time scale
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
. The topography creates variations of the thermal gradient on the liquid side which drives a convective heat transfer in the liquid acting to erase the topography by transporting the latent heat released by freezing in topography lows to topography highs where melting occurs. Other sources of motions in the liquid can also contribute to this lateral heat transfer, which happens on a time scale
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}$
, the expression for which is derived in § 2. The ratio of the two time scales,
$\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
, controls the behaviour of the boundary. For a large value of
$\unicode[STIX]{x1D6F7}$
, the topography is set by the balance between the viscous stress in the solid and the buoyancy of the topography, the phase change acting on too long a time scale to affect the classical behaviour of the free surface. The buoyancy of the topography is responsible for making the vertical velocity drop to zero at the interface, which leads to an effectively non-penetrating boundary condition. On the other hand, for low values of
$\unicode[STIX]{x1D6F7}$
, the topography is erased by freezing and melting at a rate greater than the rate at which it is generated. The removal of the associated buoyancy leads to a non-null velocity across the interface.
This situation has already been considered in the case of the dynamics of Earth’s inner core (Alboussière, Deguen & Melzani Reference Alboussière, Deguen and Melzani2010; Monnereau et al.
Reference Monnereau, Calvet, Margerin and Souriau2010; Deguen, Alboussière & Cardin Reference Deguen, Alboussière and Cardin2013; Mizzon & Monnereau Reference Mizzon and Monnereau2013), which is the solid iron sphere at the centre of the liquid iron core of the Earth. Deguen et al. (Reference Deguen, Alboussière and Cardin2013) have derived a general formulation of the boundary condition for arbitrary values of
$\unicode[STIX]{x1D6F7}$
and shown that the application of this boundary condition to a sphere considerably changes the dynamics by decreasing the critical Rayleigh number for the onset of thermal convection and allowing a new mode of convection, the translation mode, where no deformation occurs in the sphere, melting occurs at the boundary of the advancing hemisphere, and freezing occurs at the trailing boundary.
A similar situation arises for the ice shell of some satellites of giant planets in the solar system which are believed to host a liquid ocean below their ice layer (Khurana et al. Reference Khurana, Kivelson, Stevenson, Schubert, Russell, Walker and Polanskey1998; Pappalardo et al. Reference Pappalardo, Head, Greeley, Sullivan, Pilcher, Schubert, Moore, Carr, Moore, Belton and Goldsby1998; Gaidos & Nimmo Reference Gaidos and Nimmo2000; Tobie, Choblet & Sotin Reference Tobie, Choblet and Sotin2003; Soderlund et al. Reference Soderlund, Schmidt, Wicht and Blankenship2014; Čadek et al. Reference Čadek, Tobie, Van Hoolst, Massé, Choblet, Lefèvre, Mitri, Baland, Běhounková, Bourgeois and Trinh2016). Some of the largest of such satellites can also have a layer of high-pressure ices below their ocean (Grasset, Sotin & Deschamps Reference Grasset, Sotin and Deschamps2000; Sohl et al. Reference Sohl, Hussmann, Schwentker, Spohn and Lorenz2003; Baland et al. Reference Baland, Tobie, Lefèvre and van Hoolst2014). Another situation that implies such a melt–solid interface arises on all terrestrial planets in their early stage, when their silicate layer is completely or largely molten owing to the high energy of their accretion (Solomatov Reference Solomatov2007; Elkins-Tanton Reference Elkins-Tanton2012). Convection can start in the solid mantle during its crystallisation from the magma ocean, while a liquid layer persists above and/or below (Labrosse, Hernlund & Coltice Reference Labrosse, Hernlund and Coltice2007). It is therefore interesting to consider convection in a layer, not a full sphere, when a phase-change boundary condition applies at either or both its horizontal boundaries.
Deguen (Reference Deguen2013) performed such a study in the case of a spherical shell with a central gravity linearly varying with radial position and showed that, again, a translation mode is possible and favoured in the linear stability analysis if both the upper and lower boundaries allow an easy phase change – that is if each has a low value of the
$\unicode[STIX]{x1D6F7}$
parameter. The purpose of the present paper is to extend the analysis to the plane layer situation and perform the linear stability and weakly nonlinear analysis as a function of the phase-change parameters of both horizontal boundaries.
The boundary conditions are presented in § 2; § 3 presents the translation mode of convection, § 4 presents the linear and weakly nonlinear analysis in the case when both horizontal boundaries have the same value of the phase-change parameter and § 5 shows the case when phase change is allowed only on one boundary.
2 Conservation equations and boundary conditions
We consider a layer of creeping solid that behaves like a Newtonian fluid on long time scales and that is bounded above or below or both by a liquid related to the solid by a phase change (figure 1). The temperature field at rest is the solution of the thermal conduction problem with temperatures at the boundaries,
$T^{+}$
at the top and
$T^{-}$
at the bottom, that each equal the melting temperature
$T_{m}$
at the relevant pressure. Pressure, in the context of planetary interiors, is largely dominated by the hydrostatic part. The melting temperature therefore mainly depends on the vertical coordinate. The possibility of crossing the melting temperature at both the top and bottom of our computational domain requires either a nonlinear dependence of
$T_{m}$
on pressure or, more easily, a compositional difference between the solid and both upper and lower liquid layers (Labrosse et al.
Reference Labrosse, Hernlund and Coltice2007). For simplicity here, we do not consider the dynamical effects of compositional variations. The vertical dependence of the melting temperature is linearised around the reference positions of the boundaries, owing to the smallness of their topographies compared to the total thickness of the layer,
$d$
.
Definition of the topography (exaggerated here for clarity) and the temperature for the boundary conditions. The dash-dotted lines are the reference positions for the conductive motionless solutions of the top and bottom boundaries. The plot on the right shows the reference temperature profile (thick solid line) intersecting the melting temperature at the top and bottom (thin solid lines) at temperatures
$T^{+}$
and
$T^{-}$
, respectively. Lateral variations of the topography make the intersection deviate laterally in temperature. Representative temperature profiles in the liquid sides are shown as dashed lines. In the context of planetary applications, the temperature profiles should be interpreted as deviations from the isentropic reference.
The conduction temperature profile that is used as reference is written as
$$\begin{eqnarray}\displaystyle & \displaystyle T_{0}={\displaystyle \frac{T^{+}+T^{-}}{2}}+{\displaystyle \frac{z}{d}}(T^{+}-T^{-}), & \displaystyle\end{eqnarray}$$
the reference for the vertical position
$z$
being at the centre of the domain. Deviations from the conduction temperature profiles are made dimensionless using
$\unicode[STIX]{x0394}T=T^{-}-T^{+}$
as reference and denoted by
$\unicode[STIX]{x1D703}$
. In the following, superscripts
$^{+}$
and
$^{-}$
are used for quantities pertaining to the top and bottom boundaries, respectively, and omitted in equations that apply to both boundaries.
The crossing positions of the conduction solution with the melting temperature at the top and bottom are used as the reference around which a topography height,
$h^{+}$
and
$h^{-}$
, is defined for each boundary, respectively (figure 1). These topographies can have either sign, positive upward, and need not average to 0, as will be shown below. At each phase-change interface, two thermal boundary conditions are necessary to account for the moving boundary (Crank Reference Crank1984). The temperature must equal the phase-change temperature and the heat flux discontinuity across the interface must balance the release or consumption of latent heat,
$L$
(Stefan condition). The two thermal boundary conditions are written as
$$\begin{eqnarray}\displaystyle T(h)=T_{m}(h), & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}Lv_{\unicode[STIX]{x1D719}}=\unicode[STIX]{x27E6}q\unicode[STIX]{x27E7}, & & \displaystyle\end{eqnarray}$$
with
$v_{\unicode[STIX]{x1D719}}$
the freezing rate,
$\unicode[STIX]{x1D70C}_{s}$
the density of the solid and
$\unicode[STIX]{x27E6}q\unicode[STIX]{x27E7}$
the heat flux difference between the liquid and the solid sides. These boundary conditions apply to the deformed interface and need to be projected to the reference level that is used as boundary for the computation domain. Developing equation (2.2) to first order in
$h$
gives
$$\begin{eqnarray}\displaystyle & \displaystyle T\left(\pm {\displaystyle \frac{d}{2}}\right)=T^{\pm }+\left({\displaystyle \frac{\unicode[STIX]{x2202}T_{m}^{\pm }}{\unicode[STIX]{x2202}z}}-{\displaystyle \frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}z}}\right)h^{\pm }. & \displaystyle\end{eqnarray}$$
In dimensionless form, (2.4) is written as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}\left(\pm {\displaystyle \frac{1}{2}}\right)=\left(1+{\displaystyle \frac{d}{\unicode[STIX]{x0394}T}}{\displaystyle \frac{\unicode[STIX]{x2202}T_{m}^{\pm }}{\unicode[STIX]{x2202}z}}\right){\displaystyle \frac{h^{\pm }}{d}}. & \displaystyle\end{eqnarray}$$
In the following, we assume
$h^{\pm }/d$
to be small and we apply
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}=0,\quad z=\pm {\displaystyle \frac{1}{2}}. & \displaystyle\end{eqnarray}$$
Turning to the second thermal boundary condition, the discontinuity of heat flow on the right-hand side of (2.3) is assumed to be dominated by the convective heat flow on the low-viscosity liquid side,
$f\sim \unicode[STIX]{x1D70C}_{l}c_{pl}u_{l}\unicode[STIX]{x1D6FF}T_{l}$
, with
$\unicode[STIX]{x1D70C}_{l}$
and
$c_{pl}$
the density and heat capacity of the liquid,
$u_{l}$
the characteristic liquid velocity and
$\unicode[STIX]{x1D6FF}T_{l}$
the temperature difference between the boundary and the bulk of the liquid. This difference results from variations of the topography (figure 1) and the vertical gradient of the melting temperature so that
$$\begin{eqnarray}\displaystyle & \displaystyle f\sim -\unicode[STIX]{x1D70C}_{l}c_{pl}u_{l}\left|{\displaystyle \frac{\unicode[STIX]{x2202}T_{m}}{\unicode[STIX]{x2202}z}}\right|h. & \displaystyle\end{eqnarray}$$
The temperature difference
$h\unicode[STIX]{x2202}T_{m}/\unicode[STIX]{x2202}z$
is negligible on the solid side, but crucial for the convective heat flux on the liquid side. Figure 1 shows as dashed lines the typical local temperature profiles on the liquid side of each boundary for topography highs and lows, indicating that the implied lateral variations of heat flux density should lead to melting of regions where the solid protrudes in the liquid, with freezing in depressed regions, tending to erase the topography. This behaviour is ensured by the anti-correlation of
$f$
and
$h$
in (2.7), independently of the sign of
$\unicode[STIX]{x2202}T_{m}/\unicode[STIX]{x2202}z$
, and this applies to both top and bottom boundaries. The case of
$\unicode[STIX]{x2202}T_{m}/\unicode[STIX]{x2202}z<0$
depicted here for the top boundary is the most usual and the opposite case depicted here for the bottom boundary is encountered for water. Note, however, that in the context of planetary applications, the temperature considered here in the liquid layers and depicted on figure 1 is in fact the deviation from the reference isentropic temperature profile (Jeffreys Reference Jeffreys1930; Deguen et al.
Reference Deguen, Alboussière and Cardin2013) and the pressure derivative of the actual melting temperature need not be negative to have a liquid underlying the solid layer. Assuming that the convective heat flow on the liquid side dominates the right-hand side of (2.3), we write
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{s}Lv_{\unicode[STIX]{x1D719}}\sim -\unicode[STIX]{x1D70C}_{l}c_{pl}u_{l}\left|{\displaystyle \frac{\unicode[STIX]{x2202}T_{m}}{\unicode[STIX]{x2202}z}}\right|h. & \displaystyle\end{eqnarray}$$
The freezing rate is related to the vertical velocity
$w$
across the boundary and the rate of change of the topography as
$$\begin{eqnarray}\displaystyle & \displaystyle v_{\unicode[STIX]{x1D719}}^{\pm }=\pm {\displaystyle \frac{\unicode[STIX]{x2202}h^{\pm }}{\unicode[STIX]{x2202}t}}\mp w. & \displaystyle\end{eqnarray}$$
Combining with (2.8) gives
$$\begin{eqnarray}\displaystyle & \displaystyle w\mp {\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}}={\displaystyle \frac{\unicode[STIX]{x1D70C}_{l}c_{pl}u_{l}}{\unicode[STIX]{x1D70C}_{s}L}}\left|{\displaystyle \frac{\unicode[STIX]{x2202}T_{m}}{\unicode[STIX]{x2202}z}}\right|h\equiv {\displaystyle \frac{h}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}}}, & \displaystyle\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}={\displaystyle \frac{\unicode[STIX]{x1D70C}_{s}L}{\unicode[STIX]{x1D70C}_{l}c_{pl}u_{l}\left|{\displaystyle \frac{\unicode[STIX]{x2202}T_{m}}{\unicode[STIX]{x2202}z}}\right|}} & \displaystyle\end{eqnarray}$$
the characteristic phase-change time scale for changing the topography by transferring latent heat from regions where it is released to places where it is consumed.
$u_{l}$
depends on the dynamics of the liquid, which is not solved in this paper. The uncertainty in this quantity as well as the scaling coefficients implied by the
${\sim}$
sign in (2.7) and (2.8) are all combined to make
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}$
the control parameter in our study.
Across the boundaries, the total traction must be continuous. Assuming that the topography is small (i.e. the horizontal gradient of
$h^{\pm }$
is small compared to 1,
$|\unicode[STIX]{x1D6FB}_{h}h^{\pm }|\ll 1$
), for the vertical component this is written as
$$\begin{eqnarray}\displaystyle & \displaystyle -P_{s}(h^{\pm })+2\unicode[STIX]{x1D702}{\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}}=-P_{l}(h^{\pm }), & \displaystyle\end{eqnarray}$$
where
$P$
is total pressure,
$s$
and
$l$
are for the solid and liquid sides, respectively, and
$\unicode[STIX]{x1D702}$
is the dynamic viscosity of the solid. The total pressure on the solid and liquid sides is split into its hydrostatic part,
$P(0)-\unicode[STIX]{x1D70C}_{s,l}gh^{\pm }$
(
$z=0$
being the reference for
$h$
at each boundary) and the dynamic part
$p$
. On the liquid side, viscous stress and pressure fluctuations are neglected. With these assumptions, we get
$$\begin{eqnarray}\displaystyle & \displaystyle -p+(\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{l}^{\pm })gh^{\pm }+2\unicode[STIX]{x1D702}{\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}}=0. & \displaystyle\end{eqnarray}$$
Note that the density difference across the phase-change boundary,
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\pm }=\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{l}^{\pm }$
, takes different signs at the top and bottom since the solid must be denser than the overlying liquid but less dense that the underlying one. Therefore
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{+}>0$
and
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{-}<0$
.
The topography at each boundary is produced as a result of total stress in the solid, with a typical time scale
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D702}/|\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\pm }|gd$
(the post-glacial rebound time scale, Turcotte & Schubert Reference Turcotte and Schubert2001), and erased by melting and freezing, as discussed above, with a time scale
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}$
. Both time scales are generally much shorter than the time scale for convection in the whole domain, so that we assume that the topography adjusts instantaneously to the competition between viscous stress and phase change. Therefore, we neglect
$\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}t$
in (2.10) and, combining it with (2.13) to eliminate
$h^{\pm }$
, we get
$$\begin{eqnarray}\displaystyle & \displaystyle -p+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\pm }g\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}^{\pm }w+2\unicode[STIX]{x1D702}{\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}}=0. & \displaystyle\end{eqnarray}$$
Introducing the phase-change dimensionless number (Deguen Reference Deguen2013; Deguen et al. Reference Deguen, Alboussière and Cardin2013)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}^{\pm }={\displaystyle \frac{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}^{\pm }}|\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\pm }|gd}{\unicode[STIX]{x1D702}}} & \displaystyle\end{eqnarray}$$
equation (2.14) takes the dimensionless form
$$\begin{eqnarray}\displaystyle & \displaystyle \pm \unicode[STIX]{x1D6F7}^{\pm }w+2{\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}}-p=0,\quad z=\pm {\displaystyle \frac{1}{2}}. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D6F7}^{\pm }$
is the ratio of the phase-change time scale to the viscous deformation time scale. For large values of this parameter, the boundary condition (2.16) reduces to the usual non-penetration condition,
$w=0$
, while for small values it allows a non-zero mass flow through the boundary. The physical interpretation is straightforward: if
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\ll \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}^{\pm }}$
, topography evolves without the possibility of the phase change happening and is limited by its own weight, which has to be supported by viscous stress in the solid. In practice, this means that the flow velocity goes to zero at the free interface and is very small at the reference boundaries
$z=\pm 1/2$
, which is usually modelled as a non-penetrating boundary. In the other limiting case,
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\gg \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}^{\pm }}$
, topography is removed by phase change as fast as it is created by viscous stresses, and this allows a flow across the boundary.
The liquid is assumed inviscid and therefore exerts no shear stress on the convecting solid. The topography of the boundary is assumed to be small and we approximate the horizontal component of the continuity condition for traction by a free-slip boundary condition at both horizontal boundaries,
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}}+{\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}x}}=0,\quad z=\pm {\displaystyle \frac{1}{2}}. & \displaystyle\end{eqnarray}$$
The dimensionless equations for the conservation of momentum, mass and energy are written in the classical Boussinesq approximation as
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{1}{Pr}}\left({\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}t}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+\mathit{Ra}\unicode[STIX]{x1D703}\hat{\boldsymbol{z}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}=w+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
where
$Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$
is the Prandtl number, with
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D705}$
the momentum and thermal diffusivities,
$\boldsymbol{v}=(u,v,w)$
is the fluid velocity,
$p$
is the dynamic pressure,
$\mathit{Ra}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}Tgd^{3}/\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}$
is the Rayleigh number, with
$\unicode[STIX]{x1D6FC}$
the thermal expansion coefficient, and
$\hat{\boldsymbol{z}}$
is the upward vertical unit vector. These equations have been made dimensionless using the thickness of the layer
$d$
as length scale and the thermal diffusion time
$d^{2}/\unicode[STIX]{x1D705}$
as time scale.
Since we are concerned here with convection in solid, albeit creeping, layers, we will generally consider the Prandtl number to be infinite in most of the calculations below.
3 The translation mode
The boundary condition (2.16) discussed in the previous section permits a non-zero vertical velocity across the boundaries. If both boundaries are semi-permeable (finite values of both
$\unicode[STIX]{x1D6F7}^{+}$
and
$\unicode[STIX]{x1D6F7}^{-}$
), the possibility of a uniform vertical translation arises. This situation has been explored systematically in the context of the dynamics of Earth’s inner core (Alboussière et al.
Reference Alboussière, Deguen and Melzani2010; Deguen et al.
Reference Deguen, Alboussière and Cardin2013; Mizzon & Monnereau Reference Mizzon and Monnereau2013) and in spherical shells (Deguen Reference Deguen2013) but, in the case of a spherical geometry, the horizontally average vertical velocity is still null for a translation mode. Here we show that a translation mode with a uniform vertical velocity also exists in the case of a plane layer.
We search for a solution that is independent of the horizontal direction and therefore only has a vertical velocity,
$\boldsymbol{v}=w\hat{\boldsymbol{z}}$
. The mass conservation equation (2.19) implies that
$w$
is independent of
$z$
and we consider two situations, the linear stability problem for which
$w=W\text{e}^{\unicode[STIX]{x1D70E}t}$
and the steady-state case for which
$w$
is constant. Similarly, we can write the temperature as
$\unicode[STIX]{x1D703}(z,t)=\unicode[STIX]{x1D6E9}(z)\text{e}^{\unicode[STIX]{x1D70E}t}$
to study the onset of convection in that mode, with
$\unicode[STIX]{x1D703}$
as a function of
$z$
only at steady state, and use a similar convention for pressure as
$p$
and
$P$
.
3.1 Linear stability analysis
The conservation equations (2.18)–(2.20) linearised around the hydrostatic state reduce to two equations:
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D70E}}{Pr}}W=-\text{D}P+\mathit{Ra}\unicode[STIX]{x1D6E9}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}\unicode[STIX]{x1D6E9}=W+\text{D}^{2}\unicode[STIX]{x1D6E9}, & \displaystyle\end{eqnarray}$$
with
$\text{D}\equiv \text{d}/\text{d}z$
. For neutral stability,
$\unicode[STIX]{x1D70E}=0$
, solving in turn (3.2) for
$\unicode[STIX]{x1D6E9}$
and (3.1) for
$P$
subject to the boundary conditions (2.6) and (2.16) leads to
$$\begin{eqnarray}\displaystyle & \displaystyle [\mathit{Ra}-12(\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})]W=0. & \displaystyle\end{eqnarray}$$
A non-trivial solution for
$W$
can then exist for
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}=\mathit{Ra}_{c}=12(\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-}), & \displaystyle\end{eqnarray}$$
which is the condition for marginal stability of the translation mode.
This system of equations can also be solved for a finite value of
$\unicode[STIX]{x1D70E}$
in order to relate it to
$\mathit{Ra}$
. Equation (3.2) subject to boundary conditions
$\unicode[STIX]{x1D703}(\pm 1/2)=0$
gives
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}={\displaystyle \frac{W}{\unicode[STIX]{x1D70E}}}\left[1-2{\displaystyle \frac{\sinh (\unicode[STIX]{x1D70E}^{1/2}/2)}{\sinh (\unicode[STIX]{x1D70E}^{1/2})}}\cosh (\unicode[STIX]{x1D70E}^{1/2}z)\right]. & \displaystyle\end{eqnarray}$$
Inserting this expression into (3.1) and solving for
$P$
, we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle P=cst+\left({\displaystyle \frac{\mathit{Ra}}{\unicode[STIX]{x1D70E}}}-{\displaystyle \frac{\unicode[STIX]{x1D70E}}{Pr}}\right)Wz-2\mathit{Ra}\,W\unicode[STIX]{x1D70E}^{-3/2}{\displaystyle \frac{\sinh (\unicode[STIX]{x1D70E}^{1/2}/2)}{\sinh (\unicode[STIX]{x1D70E}^{1/2})}}\sinh (\unicode[STIX]{x1D70E}^{1/2}z). & \displaystyle\end{eqnarray}$$
Using the boundary condition (2.16) at
$z=1/2$
allows us to determine the integration constant, which gives
$$\begin{eqnarray}\displaystyle \displaystyle P & = & \displaystyle \unicode[STIX]{x1D6F7}^{+}W+\left({\displaystyle \frac{\mathit{Ra}}{\unicode[STIX]{x1D70E}}}-{\displaystyle \frac{\unicode[STIX]{x1D70E}}{Pr}}\right)W(z-1/2)\nonumber\\ \displaystyle & & \displaystyle -\,2\mathit{Ra}\,W\unicode[STIX]{x1D70E}^{-3/2}{\displaystyle \frac{\sinh (\unicode[STIX]{x1D70E}^{1/2}/2)}{\sinh (\unicode[STIX]{x1D70E}^{1/2})}}[\sinh (\unicode[STIX]{x1D70E}^{1/2}z)-\sinh (\unicode[STIX]{x1D70E}^{1/2}/2)].\end{eqnarray}$$
Finally, using the boundary condition at
$z=-1/2$
,
$-\unicode[STIX]{x1D719}^{-}W=P(-1/2)$
, gives, after rearranging, the following dispersion equation:
$$\begin{eqnarray}\displaystyle & \displaystyle 0={\displaystyle \frac{\unicode[STIX]{x1D70E}^{2}}{Pr(\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})}}+\unicode[STIX]{x1D70E}+{\displaystyle \frac{\mathit{Ra}}{\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-}}}\left[2\unicode[STIX]{x1D70E}^{-1/2}{\displaystyle \frac{\cosh \unicode[STIX]{x1D70E}^{1/2}-1}{\sinh \unicode[STIX]{x1D70E}^{1/2}}}-1\right]. & \displaystyle\end{eqnarray}$$
An approximate solution for small
$\unicode[STIX]{x1D70E}$
can be obtained by developing the ratio of
$\cosh$
and
$\sinh$
functions to the second order in
$\unicode[STIX]{x1D70E}$
, which gives
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}={\displaystyle \frac{10}{1+{\displaystyle \frac{120}{Pr\,\mathit{Ra}}}}}\left(1-{\displaystyle \frac{12(\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})}{\mathit{Ra}}}\right). & \displaystyle\end{eqnarray}$$
The critical Rayleigh number, obtained by setting
$\unicode[STIX]{x1D70E}=0$
, is the same as that of (3.4). If
$Gr_{T}\equiv Pr\,\mathit{Ra}$
(similar to the Grashof number but with
$\unicode[STIX]{x1D705}$
in place of
$\unicode[STIX]{x1D708}$
) is large, the expression for the growth rate reduces to
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}=10\left(1-{\displaystyle \frac{12(\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})}{\mathit{Ra}}}\right). & \displaystyle\end{eqnarray}$$
In the limit of a large
$\unicode[STIX]{x1D70E}$
,
$$\begin{eqnarray}\displaystyle & \displaystyle 2\unicode[STIX]{x1D70E}^{-1/2}{\displaystyle \frac{\cosh \unicode[STIX]{x1D70E}^{1/2}-1}{\sinh \unicode[STIX]{x1D70E}^{1/2}}}-1\rightarrow -1 & \displaystyle\end{eqnarray}$$
and the dispersion relation reduces to
$$\begin{eqnarray}\displaystyle & \displaystyle 0={\displaystyle \frac{\unicode[STIX]{x1D70E}^{2}}{Gr_{T}}}+{\displaystyle \frac{\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-}}{\mathit{Ra}}}\unicode[STIX]{x1D70E}-1. & \displaystyle\end{eqnarray}$$
The positive root is
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}={\displaystyle \frac{\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-}}{\mathit{Ra}}}{\displaystyle \frac{Gr_{T}}{2}}\left[\sqrt{1+{\displaystyle \frac{4}{Gr_{T}}}\left({\displaystyle \frac{\mathit{Ra}}{\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-}}}\right)^{2}}-1\right] & \displaystyle\end{eqnarray}$$
which reduces to
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}={\displaystyle \frac{\mathit{Ra}}{\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-}}} & \displaystyle\end{eqnarray}$$
in the limit of
$(1/Gr_{T})(\mathit{Ra}/\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})^{2}\ll 1$
. The growth rate in the large-
$Gr_{T}$
limit is plotted as function of
$\mathit{Ra}/\mathit{Ra}_{c}$
on figure 2.
Instability growth rate
$\unicode[STIX]{x1D70E}$
as a function of
$\mathit{Ra}/\mathit{Ra}_{c}$
, for infinite
$Gr_{T}$
, as given by the numerical solution of the full dispersion relation (solid blue line), and by the small- and large-
$\unicode[STIX]{x1D70E}$
approximations (black dashed lines).
3.2 Steady-state translation
The steady-state finite-amplitude translation mode is solution of
$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\text{D}p+\mathit{Ra}\unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle w\text{D}\unicode[STIX]{x1D703}=w+\text{D}^{2}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$
Solving first the energy balance equation (3.16) subject to boundary conditions (2.6) gives
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}=z+{\displaystyle \frac{\cosh \left({\displaystyle \frac{w}{2}}\right)-\text{e}^{wz}}{2\sinh \left({\displaystyle \frac{w}{2}}\right)}}\Rightarrow T={\displaystyle \frac{1}{2}}+{\displaystyle \frac{\cosh \left({\displaystyle \frac{w}{2}}\right)-\text{e}^{wz}}{2\sinh \left({\displaystyle \frac{w}{2}}\right)}}. & \displaystyle\end{eqnarray}$$
Using the momentum balance equation (3.15) and the boundary conditions (2.16) then gives
$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})w=\mathit{Ra}\left[{\displaystyle \frac{\cosh \left({\displaystyle \frac{w}{2}}\right)}{2\sinh \left({\displaystyle \frac{w}{2}}\right)}}-{\displaystyle \frac{1}{w}}\right]. & \displaystyle\end{eqnarray}$$
This transcendental equation relates the translation velocity
$w$
to the Rayleigh number.
Close to onset, assuming the Péclet number,
$|w|$
, to be small, equation (3.18) can be developed as a function of
$(\mathit{Ra}-\mathit{Ra}_{c})/\mathit{Ra}_{c}$
to give to leading order
$$\begin{eqnarray}\displaystyle & \displaystyle w=\pm 2\sqrt{15{\displaystyle \frac{\mathit{Ra }-\mathit{Ra }_{c}}{\mathit{Ra}_{c}}}}. & \displaystyle\end{eqnarray}$$
The corresponding temperature anomaly is
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}={\displaystyle \frac{w}{8}}(1-4z^{2})+O(w^{2}), & \displaystyle\end{eqnarray}$$
showing that the temperature only differs from the conduction solution by an amount proportional to the Péclet number.
For a large Péclet number,
$|w|\gg 1$
, equation (3.18) reduces to
$$\begin{eqnarray}\displaystyle & \displaystyle w\sim \pm {\displaystyle \frac{\mathit{Ra}}{2(\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})}}=\pm {\displaystyle \frac{6\mathit{Ra}}{\mathit{Ra}_{c}}}. & \displaystyle\end{eqnarray}$$
Figure 3 shows how the translation velocity
$|w|$
depends on Rayleigh number, computed using the full equation (3.18) and either the low- or the large-velocity development. It shows that the transition between the two regimes happens for
$\mathit{Ra}\sim 2\mathit{Ra}_{c}$
.
In the high-Péclet-number regime, the temperature anomaly takes a simple form:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}\sim z+\text{sgn}(w)\left[{\displaystyle \frac{1}{2}}-\text{e}^{w(z-\text{sgn}(w)/2)}\right]\Rightarrow T\sim {\displaystyle \frac{1}{2}}[1+\text{sgn}(w)]-\text{sgn}(w)\text{e}^{w(z-\text{sgn}(w)/2)}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The exponential in the last equation is negligible everywhere except close to the upper boundary (
$z=1/2$
; respectively lower boundary,
$z=-1/2$
) when
$w\gg 1$
(respectively
$w\ll -1$
). Therefore, the temperature is essentially equal to that imposed at the boundary the fluid originates from (
$0$
at the top,
$1$
at the bottom) and adjusts to that of the opposite side in a boundary layer of thickness
$\unicode[STIX]{x1D6FF}\sim 1/w$
. In dimensional units,
$\unicode[STIX]{x1D6FF}$
is simply defined as the thickness that makes the Péclet number approximately 1:
$\mathit{Pe}=w\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D705}\sim 1$
. Figure 4 shows the temperature profiles for the upward and downward translation modes computed both with the exact (3.17) and approximate (3.22) expressions, showing that the approximation is quite good.
Temperature profile in the translation mode for
$(\mathit{Ra}-\mathit{Ra}_{c})/\mathit{Ra}_{c}=5$
. The solid (respectively dashed) line is for the ascending (respectively descending) mode calculated using the full equation (3.17) and the up (respectively down) triangles are obtained using the approximate equation (3.22).
The steady-state velocity given by (3.21) can also be obtained from a simple physical argument. In the steady-translation regime, the (uniform) topography at each boundary is related to the translation velocity and the phase-change time scale by
$$\begin{eqnarray}\displaystyle & \displaystyle h^{\pm }=\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}^{\pm }}w. & \displaystyle\end{eqnarray}$$
In steady state, the excess (respectively deficit) weight of the cooler (respectively warmer) solid layer is balanced by the sum of pressure deviations from the hydrostatic equilibrium at both boundaries as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{0}g{\displaystyle \frac{\unicode[STIX]{x0394}Td}{2}}=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{+}gh^{+}+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{-}gh^{-}, & \displaystyle\end{eqnarray}$$
where the temperature in the solid layer has been assumed uniform, i.e. the contribution of the boundary layer to its buoyancy has been neglected. This gives for the translation velocity
$$\begin{eqnarray}\displaystyle & \displaystyle w={\displaystyle \frac{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{0}g\unicode[STIX]{x0394}Td}{2(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{+}g\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}^{+}}+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{-}g\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}^{-}})}}. & \displaystyle\end{eqnarray}$$
In dimensionless form, this is exactly (3.21).
It is also worth considering the heat transfer efficiency in the translation mode. Equation (3.16) can be integrated to show that
$wT-\text{D}T$
is independent of
$z$
and this implies that
$w=\text{D}T(-1/2)-\text{D}T(1/2)$
, meaning that the difference between the conductive heat fluxes across the horizontal boundaries is equal to the advection by translation. Figure 4 shows that the heat flow (Nusselt number
$\mathit{Nu}$
) should be computed on the exit side, where a boundary layer is produced:
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Nu}=-\text{D}T\left(\text{sgn}(w){\displaystyle \frac{1}{2}}\right)=|w|-\text{D}T\left(-\text{sgn}(w){\displaystyle \frac{1}{2}}\right)=|w|+{\displaystyle \frac{w\text{e}^{-|w|/2}}{2\sinh (w/2)}}. & \displaystyle\end{eqnarray}$$
The small- and large-
$|w|$
limit cases give
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Nu}=1+{\displaystyle \frac{|w|}{2}}=1+\sqrt{15{\displaystyle \frac{\mathit{Ra }-\mathit{Ra }_{c}}{\mathit{Ra}_{c}}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Nu}=|w|=6{\displaystyle \frac{\mathit{Ra}}{\mathit{Ra }_{c}}}, & \displaystyle\end{eqnarray}$$
respectively. The large-Rayleigh-number behaviour is in striking contrast to the situation encountered for standard Rayleigh–Bénard convection, for which
$\mathit{Nu}\sim \mathit{Ra}^{\unicode[STIX]{x1D6FD}}$
with
$\unicode[STIX]{x1D6FD}\sim 1/3$
.
4 Non-translating modes with
$\unicode[STIX]{x1D6F7}^{+}=\unicode[STIX]{x1D6F7}^{-}$
In this section, we consider the situation with values of the phase-change parameter of both boundaries equal,
$\unicode[STIX]{x1D6F7}\equiv \unicode[STIX]{x1D6F7}^{+}=\unicode[STIX]{x1D6F7}^{-}$
.
4.1 Linear stability
Non-translating solutions can be obtained using standard approaches for the classical Rayleigh–Bénard problem. For the linear stability problem, a solution using separation of variables is sought, i.e.
$u=U(z)\text{e}^{\text{i}kx}\text{e}^{\unicode[STIX]{x1D70E}t}$
, and similarly for
$w$
,
$p$
and
$\unicode[STIX]{x1D703}$
. Linearised equations (2.18) to (2.20) reduce to
$$\begin{eqnarray}\displaystyle & \displaystyle \text{i}kU+\text{D}W=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Pr[-\text{i}kP+(\text{D}^{2}-k^{2})U]=\unicode[STIX]{x1D70E}U, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Pr[-\text{D}P+(\text{D}^{2}-k^{2})W+\mathit{Ra}\unicode[STIX]{x1D6E9}]=\unicode[STIX]{x1D70E}W, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle W+(\text{D}^{2}-k^{2})\unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6E9} & \displaystyle\end{eqnarray}$$
since, at the linear stage, the problem is fully degenerate in terms of orientation of the mode, which can be taken as depending only on
$x$
. These equations must be complemented by boundary conditions applying at
$z=\pm 1/2$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \text{D}U+\text{i}kW=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \pm \unicode[STIX]{x1D6F7}^{\pm }W+2\text{D}W-P=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}=0. & \displaystyle\end{eqnarray}$$
This forms a generalised eigenvalue problem that we solve using a Chebyshev-collocation pseudo-spectral approach (e.g. Canuto et al.
Reference Canuto, Hussaini, Quarteroni and Zang1988; Guo, Labrosse & Narayanan Reference Guo, Labrosse and Narayanan2012). Given the Chebyshev–Gauss–Lobatto nodal point
$z_{i}=\cos (i\unicode[STIX]{x03C0}/N),i=0\ldots N$
, in the interval
$[-1,1]$
, the values of the
$z$
-dependent mode functions at
$z_{i}/2$
are denoted as
$U_{i}$
for
$U$
and similarly for other variables. Division by two is required here to map the interval on which Chebyshev polynomials are defined onto
$[-1/2,1/2]$
. The
$k$
th derivative of each function at the nodal points is related to the nodal values of the function itself by differentiation matrices:
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{U}^{(k)}=\unicode[STIX]{x1D63F}^{(k)}\boldsymbol{\cdot }\boldsymbol{U}. & \displaystyle\end{eqnarray}$$
The calculation of the differentiation matrices is done using a Python adaptation (available at https://github.com/labrosse/dmsuite) of DMSUITE (Weideman & Reddy Reference Weideman and Reddy2000). With these differentiation matrices, the system of equations (4.1) to (4.4) can be written as a generalised eigenvalue problem of the form
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}\boldsymbol{\cdot }\boldsymbol{X}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D64D}\boldsymbol{\cdot }\boldsymbol{X}, & \displaystyle\end{eqnarray}$$
with
$\boldsymbol{X}=(\boldsymbol{P};\boldsymbol{U};\boldsymbol{W};\unicode[STIX]{x1D723})^{\text{T}}$
the global vertical mode vector composed of the concatenation of vectors
$\boldsymbol{P}$
,
$\boldsymbol{U}$
,
$\boldsymbol{W}$
and
$\unicode[STIX]{x1D723}$
, and
$\unicode[STIX]{x1D647}$
and
$\unicode[STIX]{x1D64D}$
two matrices representing the system with its boundary conditions. The general structure of
$\unicode[STIX]{x1D647}$
reads as
with
$\unicode[STIX]{x1D644}$
and
$\mathbf{0}$
the identity and zero matrices, respectively. The restrictions of line and column indices, indicated on the right and top of the matrix respectively, are necessary to leave out the boundary points from applications of (4.1) to (4.4), since these follow equations (4.5)–(4.7) instead. For example, in the second line of the matrix that represents (4.5), only the first line (index
$0$
) of the matrice
$\mathbf{0},\unicode[STIX]{x1D63F},\text{i}k\unicode[STIX]{x1D644}$
and
$\mathbf{0}$
are present. Note that the boundary values for the temperature are simply left out since the Dirichlet boundary condition (4.7) is, in a collocation approach, naturally enforced by removing the extreme Chebyshev points.
The
$\unicode[STIX]{x1D64D}$
matrix contains ones on the diagonal corresponding to the interior points of the equations for
$\boldsymbol{U}$
,
$\boldsymbol{W}$
and
$\unicode[STIX]{x1D723}$
, with zeros elsewhere. When solving for an infinite Prandtl number, which is the case below, the interior points for the
$\boldsymbol{U}$
and
$\boldsymbol{W}$
equations are also set to 0, leaving ones only for the interior points of the
$\unicode[STIX]{x1D723}$
equation. The resulting system is singular and many eigenvalues are infinite, one for each zero on the diagonal of the
$\unicode[STIX]{x1D64D}$
matrix. Filtering these spurious eigenvalues leaves us with the relevant eigenvalues that are used to assess stability. For any values of
$\unicode[STIX]{x1D6F7}^{-}$
,
$\unicode[STIX]{x1D6F7}^{+}$
and
$k$
, the minimum value of
$\mathit{Ra}$
that makes the real part of one of the eigenvalues become positive is the critical Rayleigh number for perturbations with that wavenumber. Minimising
$\mathit{Ra}$
as function of
$k$
gives the critical Rayleigh number for all infinitesimal perturbations. Figure 5 shows the evolution of the critical Rayleigh number and the associated wavenumber as functions of the value of
$\unicode[STIX]{x1D6F7}^{\pm }$
, both taken equal,
$\unicode[STIX]{x1D6F7}^{+}=\unicode[STIX]{x1D6F7}^{-}=\unicode[STIX]{x1D6F7}$
. One can see that the classical value derived by Rayleigh (Reference Rayleigh1916) is recovered when
$\unicode[STIX]{x1D6F7}\rightarrow \infty$
, as expected. In the other limit,
$\unicode[STIX]{x1D6F7}\rightarrow 0$
, the critical Rayleigh number follows the analytical expression obtained for the translation mode (§ 3), while
$k\rightarrow 0$
, as expected.
Critical Rayleigh number (a) and wavenumber (b) as functions of the phase-change numbers, both taken equal here. Filled circles are results of the calculation using the Chebyshev-collocation technique, the dash-dotted lines represent the classical
$\unicode[STIX]{x1D6F7}\rightarrow \infty$
limit, the dashed line in (a) represents the result for the translating mode (3.4) and the solid lines represent the small-
$\unicode[STIX]{x1D6F7}$
leading-order development.
The behaviour of the system in the limit of small
$\unicode[STIX]{x1D6F7}$
can be obtained using a polynomial expansion of all the functions, both in
$z$
and
$\unicode[STIX]{x1D6F7}$
. Specifically, considering the symmetry of the problem around
$z=0$
, we write the temperature as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}=\mathop{\sum }_{n=0}^{N}a_{n}z^{2n}. & \displaystyle\end{eqnarray}$$
The Hermitian character of the linear problem (see appendix A) ensures that
$\unicode[STIX]{x1D70E}$
is real and, therefore,
$\unicode[STIX]{x1D70E}=0$
at onset. Then
$W$
and
$U$
can be obtained using (4.4) and (4.1). Equations (4.2) and (4.3) then provide two expressions for D
$P$
and their equality implies several equations, one for each polynomial order considered. All the functions are developed to the same order as the temperature,
$2N$
. Note that even if the definition of
$\unicode[STIX]{x1D6E9}$
for a given
$N$
only requires
$N+1$
coefficients
$a_{n}$
, the development of the other profiles to the same order requires the inclusion of
$a_{n}$
for values up to
$n=N+2$
because of the derivatives in the linear system. Using, for example,
$N=2$
gives a pressure gradient D
$P$
that contains terms in
$z^{2n},n=0\ldots 2$
, and provides therefore three independent equations for the equality between the two expressions. With the symmetry considered here, the boundary conditions (4.5)–(4.7) introduce three additional equations for the coefficients
$a_{n}$
.
Setting first
$\unicode[STIX]{x1D6F7}=0$
leads to a non-trivial solution only for
$\mathit{Ra}=0$
and
$k=0$
, the solution being equal to the low-
$\unicode[STIX]{x1D6F7}$
development of the translation solution. To go beyond that, each coefficient
$a_{n}$
is itself developed as a polynomial of
$\unicode[STIX]{x1D6F7}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle a_{n}=\mathop{\sum }_{j=0}^{J}a_{n,j}\unicode[STIX]{x1D6F7}^{j}. & \displaystyle\end{eqnarray}$$
Similarly, the critical Rayleigh number
$\mathit{Ra}_{c}$
and the square of the critical wavenumber
$k^{2}$
are developed in powers of
$\unicode[STIX]{x1D6F7}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{c}=\mathop{\sum }_{j=0}^{J}r_{j}\unicode[STIX]{x1D6F7}^{j},\quad k_{c}^{2}=\mathop{\sum }_{j=0}^{J}K_{j}\unicode[STIX]{x1D6F7}^{j}. & \displaystyle\end{eqnarray}$$
The three boundary conditions and the equations implied by the equality of the two pressure expressions are then written and solved for increasing degrees in the development in
$\unicode[STIX]{x1D6F7}$
. In practice, we restrict ourselves to
$N=J=2$
. At order 0 in
$\unicode[STIX]{x1D6F7}$
, the set of linear equations can admit a non-trivial solution only if the determinant of the implied matrix is zero, which provides two possible values of
$r_{0}$
. The lowest one admits a minimum,
$r_{0}=0$
, for
$K_{0}=0$
. This implies
$a_{2,0}=a_{3,0}=a_{4,0}=0$
and
$a_{1,0}=-4a_{0,0}$
. At order 1 in
$\unicode[STIX]{x1D6F7}$
, we get directly that
$a_{2,1}=a_{3,1}=a_{4,1}=0$
,
$a_{1,1}=-4a_{0,1}$
and
$r_{1}=24$
with no information on
$K_{1}$
. This is, however, obtained at the next order, where we find that
$K_{1}=9/32$
minimises
$r_{2}$
, which is then
$r_{2}=-81/256$
. The order-2 coefficients are also obtained as a function of
$a_{0,0}$
, which is the value of the maximum of
$\unicode[STIX]{x1D6E9}$
. These can then be used to determine the shape of the different functions
$\unicode[STIX]{x1D6E9}$
,
$W$
,
$U$
and
$P$
for small values of
$\unicode[STIX]{x1D6F7}$
. To leading order in
$\unicode[STIX]{x1D6F7}$
we get
$$\begin{eqnarray}\displaystyle & \displaystyle k_{c}={\displaystyle \frac{3}{4\sqrt{2}}}\sqrt{\unicode[STIX]{x1D6F7}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{c}=24\unicode[STIX]{x1D6F7}-{\displaystyle \frac{81}{256}}\unicode[STIX]{x1D6F7}^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}=(1-4z^{2})\unicode[STIX]{x1D6E9}_{max}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle W=8\unicode[STIX]{x1D6E9}_{max}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle U=-3\text{i}\sqrt{2\unicode[STIX]{x1D6F7}}z\unicode[STIX]{x1D6E9}_{max}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle P={\displaystyle \frac{z}{2}}(39-64z^{2})\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6E9}_{max}. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D6E9}_{max}=a_{0,0}$
is used to normalise all profiles. Note that the shape of the temperature (4.16) and vertical velocity (4.17) profiles are of order 0 in
$\unicode[STIX]{x1D6F7}$
and are equal to their counterpart in the steady-state-translation solution (3.20). The small-
$\unicode[STIX]{x1D6F7}$
development of the solution to the linear problem can be compared to the results obtained using the Chebyshev-collocation method for cross-validation. The match between the mode profiles is very good for
$\unicode[STIX]{x1D6F7}\leqslant 0.1$
. Figure 5 shows the variation of
$\mathit{Ra}_{c}$
and
$k_{c}$
as functions of
$\unicode[STIX]{x1D6F7}$
as computed by the Chebyshev-collocation approach (in solid symbols) as well as the analytical value classically obtained for non-penetrating conditions and the small-
$\unicode[STIX]{x1D6F7}$
expansion. Additionally, figure 6 shows the variation of the maximum of profiles of
$P$
,
$U$
and
$W$
, that of
$\unicode[STIX]{x1D6E9}$
being set to 1, as well as the difference between the critical Rayleigh number for uniform translation (
$24\unicode[STIX]{x1D6F7}$
) and that for a deforming mode, each as functions of
$\unicode[STIX]{x1D6F7}$
. It shows the consistency between the calculations using the Chebyshev-collocation approach and the low-
$\unicode[STIX]{x1D6F7}$
development.
Variation of the maxima of profiles of
$P$
,
$U$
and
$W$
of the first unstable mode (a–c, respectively), that for
$\unicode[STIX]{x1D6E9}$
being set to 1, as a function of
$\unicode[STIX]{x1D6F7}$
. Panel (d) shows the difference between
$24\unicode[STIX]{x1D6F7}$
and the critical Rayleigh number. On each plot, the solid circles are the results of the calculation using the Chebyshev-collocation method while the dashed lines are the low-
$\unicode[STIX]{x1D6F7}$
predictions of equations (4.15) to (4.19).
At low
$\unicode[STIX]{x1D6F7}$
, the wavelength of the first unstable mode tends to infinity as
${\sim}1/\sqrt{\unicode[STIX]{x1D6F7}}$
, which means that deformation of the solid becomes negligible. Accordingly, the viscous stress ceases to be a limiting factor for the flow, and
$Ra_{c}/\unicode[STIX]{x1D6F7}$
, which contains no viscosity, tends to a constant value. This ratio,
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\mathit{Ra}}{\unicode[STIX]{x1D6F7}}}={\displaystyle \frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}Td^{2}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\pm }\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}}}\equiv {\displaystyle \frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}_{T}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\pm }}}{\displaystyle \frac{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}}}, & \displaystyle\end{eqnarray}$$
is the ratio of the driving thermal density difference
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}_{T}$
to that involved in the phase change, times the ratio of the thermal time scale to the phase-change time scale, and can be considered as the effective Rayleigh number in the low-
$\unicode[STIX]{x1D6F7}$
limit.
First unstable mode for three different values of
$\unicode[STIX]{x1D6F7}^{+}=\unicode[STIX]{x1D6F7}^{-}$
:
$10^{5}$
(a),
$10$
(b) and
$10^{-2}$
(c). The colour represents temperature and the flow lines thickness is proportional to the norm of the velocity. Note that the
$x$
range is different in (a), (b) and (c) due to the change in wavelength. The
$z$
axis has been scaled accordingly in (a) and (b), but not in (c) as the height of the figure would be too small to read.
Figure 7 shows the first unstable mode for different values of the phase-change parameter. In the case of
$\unicode[STIX]{x1D6F7}=10^{5}$
, the critical Rayleigh number and wavenumber are very close to that obtained using classical non-penetrating boundary conditions (figure 5), and so is the first unstable mode. For
$\unicode[STIX]{x1D6F7}=10$
, the critical Rayleigh number has already decreased significantly (
$\mathit{Ra}_{c}=190$
), the critical wavelength significantly increased (
$\unicode[STIX]{x1D706}_{c}=4.55$
) and the critical mode displays streamlines that cross both boundaries. For
$\unicode[STIX]{x1D6F7}=10^{-2}$
, the critical Rayleigh number is slightly less than
$0.24$
, the critical wavelength is approximately
$115$
and streamlines are essentially vertical. At each horizontal position, this mode of convection has exactly the same shape as the linearly unstable translation mode, but it is modulated laterally, with a very long wavelength that increases as
${\sim}1/\sqrt{\unicode[STIX]{x1D6F7}}$
when
$\unicode[STIX]{x1D6F7}\rightarrow 0$
. The fact that this makes the critical Rayleigh number smaller than that for pure solid-body translation is rather mysterious.
The critical Rayleigh number for the instability for the non-null
$k$
mode is always lower than that for pure translation, as shown by (4.15) and figure 5, and should therefore always be favoured. This might be true in an infinite layer but, in practical cases, the horizontal direction is periodic, either in numerical models or in a planetary mantle. In that case, the minimal value of
$k$
that can be attained is
$2\unicode[STIX]{x03C0}/L$
, with
$L$
the horizontal periodicity. If the value of
$k$
corresponding to the critical Rayleigh number is smaller than
$2\unicode[STIX]{x03C0}/L$
, the translation mode could still be favoured. The study of the stability of the uniformly translating solution with respect to laterally varying modes is a simple extension to the stability of the conduction solution. Considering now that
$(p,\boldsymbol{v},\unicode[STIX]{x1D703})$
are infinitesimal perturbations with respect to the steady-translation solution
$(p_{t},w_{t}\hat{\boldsymbol{z}},T_{t})$
, the only equation to be modified compared to that treated in § 4.1 at infinite Prandtl number is the temperature equation, which now reads
$$\begin{eqnarray}\displaystyle & \displaystyle (\text{D}^{2}-k^{2})\unicode[STIX]{x1D6E9}-w_{t}\text{D}\unicode[STIX]{x1D6E9}-W\text{D}T_{t}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6E9} & \displaystyle\end{eqnarray}$$
instead of (4.4). Using the steady-translation solution provided in § 3.2, this equation can be implemented in the stability calculation to compute the growth rate of a deforming perturbation of wavenumber
$k$
when a steady-translation solution is in place for a given Rayleigh number above the critical value for the translation solution. We denote by
$\unicode[STIX]{x1D700}=(\mathit{Ra}-\mathit{Ra}_{c})/\mathit{Ra}_{c}$
the reduced Rayleigh number,
$\mathit{Ra}_{c}=12(\unicode[STIX]{x1D6F7}^{+}+\unicode[STIX]{x1D6F7}^{-})$
, being here the critical value for the onset of uniform translation. When
$\unicode[STIX]{x1D700}$
tends to zero, the translation velocity
$w_{t}$
tends to zero and the system of equations tends to that solved for the stability of the steady-conduction solution. But since
$\unicode[STIX]{x1D700}=0$
corresponds to the critical Rayleigh number for the translation solution, which is finitely greater than the critical value for the instability with finite
$k$
, we expect a finite instability growth rate in a finite band of wavenumbers. We therefore expect an infinitely slow translation solution to be unstable with respect to deforming modes. However, when the Rayleigh number is increased above the critical value for the translation mode, we expect this translation mode with a finite velocity to become more stable since perturbations with a finite
$k$
are then transported away by translation. Figure 8 indeed shows that, for a given value of the phase-change number
$\unicode[STIX]{x1D6F7}$
(equal for both boundaries here), increasing the Rayleigh number above the critical value for the translation mode, and therefore the steady-state-translation velocity, the linear growth rate of the deforming mode decreases. For a given Rayleigh number, the growth rate curve as a function of wavenumber displays a maximum; this maximum decreases with Rayleigh number and eventually becomes negative. There is therefore a maximum Rayleigh number beyond which the translation solution is linearly stable against any deforming perturbation. Figure 9 shows the range of unstable modes in the
$k$
–
$\unicode[STIX]{x1D700}$
space for three different values of the phase-change number. The range of Rayleigh numbers above the critical one for translation that allows the finite-
$k$
instabilities to develop shrinks when
$\unicode[STIX]{x1D6F7}$
decreases and the translation mode becomes increasingly more relevant. Figure 10 shows that the maximum growth rate of the instability at
$\unicode[STIX]{x1D700}=0$
varies linearly with
$\unicode[STIX]{x1D6F7}$
, as also does the maximum value of
$\unicode[STIX]{x1D700}$
for an instability to develop. The wavenumber for the instability is found to be equal to that for the instability of the conductive solution (figure 9), and therefore varies as
$\sqrt{\unicode[STIX]{x1D6F7}}$
(figure 5).
Growth rate of deforming perturbation over a steady-translating solution as a function of the perturbation wavenumber
$k$
, for different values of the reduced Rayleigh number
$\unicode[STIX]{x1D700}=(\mathit{Ra}-\mathit{Ra}_{c})/\mathit{Ra}_{c}$
and for
$\unicode[STIX]{x1D6F7}^{+}=\unicode[STIX]{x1D6F7}^{-}=1$
.
Range of wavenumbers as a function of the reduced Rayleigh number for which the translation solution is unstable versus deforming modes. Three different shaded regions for three different values of
$\unicode[STIX]{x1D6F7}$
are represented. For each shaded area, the dashed line represents the values of the wavenumber giving the maximum growth rate as a function of the reduced Rayleigh number.
Maximum growth rate for a non-null
$k$
mode at the critical Rayleigh number for the onset of the translation mode (a) and maximum reduced Rayleigh number for a positive growth rate of a deforming instability over a finite-amplitude translation mode (b), as functions of the phase-change number.
4.2 Weakly nonlinear analysis
Going beyond the linear stability is necessary to assess the behaviour of the system at Rayleigh numbers larger than the critical value, in particular to investigate the heat transfer efficiency of the convective system. We here follow the approach classically developed for weakly nonlinear dynamics (Malkus & Veronis Reference Malkus and Veronis1958; Schlüter, Lortz & Busse Reference Schlüter, Lortz and Busse1965; Manneville Reference Manneville2004). The system of partial differential equations (2.18)–(2.20) is separated into its linear and nonlinear parts as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}(\unicode[STIX]{x2202}_{t},\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{z},\mathit{Ra})\boldsymbol{X}=\boldsymbol{N}(\boldsymbol{X},\boldsymbol{X}), & \displaystyle\end{eqnarray}$$
with
$\boldsymbol{X}=(p;u;w;\unicode[STIX]{x1D703})^{\text{T}}$
and for an infinite Prandtl number case
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}=\left(\begin{array}{@{}cccc@{}}0 & \unicode[STIX]{x2202}_{x} & \unicode[STIX]{x2202}_{z} & 0\\ -\unicode[STIX]{x2202}_{x} & \unicode[STIX]{x1D6FB}^{2} & 0 & 0\\ -\unicode[STIX]{x2202}_{z} & 0 & \unicode[STIX]{x1D6FB}^{2} & \mathit{Ra}\\ 0 & 0 & 1 & \unicode[STIX]{x1D6FB}^{2}-\unicode[STIX]{x2202}_{t}\end{array}\right),\quad \boldsymbol{N}(\boldsymbol{X}_{l},\boldsymbol{X}_{m})=\left[\begin{array}{@{}c@{}}0\\ 0\\ 0\\ u_{l}\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D703}_{m}+w_{l}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D703}_{m}\end{array}\right]. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The linear operator is further developed around the critical Rayleigh number as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}=\unicode[STIX]{x1D647}_{c}-(\mathit{Ra}-\mathit{Ra}_{c})\unicode[STIX]{x1D648}. & \displaystyle\end{eqnarray}$$
By giving
$\mathit{Ra}_{c}$
as weight to the
$\unicode[STIX]{x1D703}$
part in the dot product
$\langle \bullet |\bullet \rangle$
, it can be shown that the operator
$\unicode[STIX]{x1D647}_{c}$
is self-adjoint (Hermitian),
$\langle \boldsymbol{X}_{2}|\unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{1}\rangle =\langle \unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{2}|\boldsymbol{X}_{1}\rangle$
(see appendix for details). Among other things, it implies that all its eigenvalues are real and the marginal state is characterised by
$\unicode[STIX]{x2202}_{t}=0$
. The solution
$\boldsymbol{X}$
and the Rayleigh number are developed as
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{X}=\unicode[STIX]{x1D716}\boldsymbol{X}_{1}+\unicode[STIX]{x1D716}^{2}\boldsymbol{X}_{2}+\unicode[STIX]{x1D716}^{3}\boldsymbol{X}_{3}+\cdots \,, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}=\mathit{Ra}_{c}+\unicode[STIX]{x1D716}\mathit{Ra}_{1}+\unicode[STIX]{x1D716}^{2}\mathit{Ra}_{2}+\cdots \,, & \displaystyle\end{eqnarray}$$
and (4.22) leads to a set of equations for increasing orders of
$\unicode[STIX]{x1D716}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{1}=\mathbf{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{2}=\boldsymbol{N}(\boldsymbol{X}_{1},\boldsymbol{X}_{1})+\mathit{Ra}_{1}\unicode[STIX]{x1D648}\boldsymbol{X}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{3}=\boldsymbol{N}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})+\boldsymbol{N}(\boldsymbol{X}_{2},\boldsymbol{X}_{1})+\mathit{Ra}_{1}\unicode[STIX]{x1D648}\boldsymbol{X}_{2}+\mathit{Ra}_{2}\unicode[STIX]{x1D648}\boldsymbol{X}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{n}=\mathop{\sum }_{l=1}^{n-1}\boldsymbol{N}(\boldsymbol{X}_{l},\boldsymbol{X}_{n-l})+\mathop{\sum }_{l=1}^{n-1}\mathit{Ra}_{l}\unicode[STIX]{x1D648}\boldsymbol{X}_{n-l}. & \displaystyle\end{eqnarray}$$
Equation (4.27) is simply that of the linear stability problem and its solution is
$\boldsymbol{X}_{1}=\boldsymbol{X}_{c}$
, which can be suitably normalised such that the maximum value of
$W$
is 1. Taking the scalar product of equations of subsequent orders by
$\boldsymbol{X}_{1}$
and making use of the Hermitian properties of
$\unicode[STIX]{x1D647}_{c}$
provides solvability conditions (Fredholm alternative) that determine the values of
$\mathit{Ra}_{i}$
. For
$\mathit{Ra}_{1}$
one gets:
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{1}=-{\displaystyle \frac{\langle \boldsymbol{X}_{1}|\boldsymbol{N}(\boldsymbol{X}_{1},\boldsymbol{X}_{1})\rangle }{\langle \boldsymbol{X}_{1}|\unicode[STIX]{x1D648}\boldsymbol{X}_{1}\rangle }}. & \displaystyle\end{eqnarray}$$
The
$x$
dependence of
$\boldsymbol{X}_{1}$
is of the form
$\text{e}^{\text{i}k_{c}x}$
, i.e.
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{X}_{1}=\boldsymbol{Z}_{1,1}(z)\text{e}^{\text{i}k_{c}x}+\text{c.c.}, & \displaystyle\end{eqnarray}$$
with
$\boldsymbol{Z}_{1,1}(z)=(\boldsymbol{P}_{1,1}(z);\boldsymbol{U}_{1,1}(z);\boldsymbol{W}_{1,1}(z);\unicode[STIX]{x1D723}_{1,1}(z))^{\text{T}}$
the vector composed of the four vertical modes for all four variables, at degree 1 of weakly nonlinear development (first index) and for the first mode in the horizontal direction (second index).
Then,
$\boldsymbol{N}(\boldsymbol{X}_{1},\boldsymbol{X}_{1})$
contains two contributions to its
$x$
dependence, one constant and one in
$\text{e}^{\text{i}2k_{c}x}$
. It is therefore orthogonal to
$\boldsymbol{X}_{1}$
and it can then be concluded that
$\mathit{Ra}_{1}=0$
. The general solution to (4.28) is the sum of the solution to the homogeneous equation and a particular solution of the equation with a right-hand side. Since we are seeking a solution
$\boldsymbol{X}_{2}$
which adds to
$\boldsymbol{X}_{1}$
, i.e. orthogonal to it, and since
$\boldsymbol{X}_{1}$
is the general solution to the homogeneous equation, only the particular solution is of interest. The
$x$
dependence of
$\boldsymbol{X}_{2}$
will contain a constant value of the form
$\boldsymbol{Z}_{2,0}(z)$
and a term of the form
$\boldsymbol{Z}_{2,2}(z)\text{e}^{\text{i}2k_{c}x}$
. Computing the scalar product of (4.29) by
$\boldsymbol{X}_{1}$
gives the value of
$\mathit{Ra}_{2}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{2}=-{\displaystyle \frac{\langle \boldsymbol{X}_{1}|\boldsymbol{N}(\boldsymbol{X}_{2},\boldsymbol{X}_{1})\rangle +\langle \boldsymbol{X}_{1}|\boldsymbol{N}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})\rangle }{\langle \boldsymbol{X}_{1}|\unicode[STIX]{x1D648}\boldsymbol{X}_{1}\rangle }}. & \displaystyle\end{eqnarray}$$
$\boldsymbol{X}_{\mathbf{2}}$
containing a term proportional to
$\text{e}^{\text{i}2k_{c}x}$
and a term independent of
$x$
,
$\boldsymbol{N}(\boldsymbol{X}_{2},\boldsymbol{X}_{1})$
and
$\boldsymbol{N}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})$
have contributions of the form
$\text{e}^{\pm \text{i}k_{c}x}$
which can resonate with
$\boldsymbol{X}_{1}$
and make
$\mathit{Ra}_{2}$
non-null. In that case, the amplitude parameter is, to leading order,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}=\sqrt{{\displaystyle \frac{\mathit{Ra }-\mathit{Ra }_{c}}{\mathit{Ra}_{2}}}}. & \displaystyle\end{eqnarray}$$
The procedure can be extended to any higher order and the general behaviour can be predicted by recursive reasoning. In particular, it is easy to show that solutions of even and odd order contain contributions to their
$x$
dependence as even and odd powers of
$\text{e}^{\text{i}k_{c}x}$
up to their order value, i.e.
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{X}_{2n}=\mathop{\sum }_{l=0}^{n}\boldsymbol{Z}_{2n,2l}(z)\text{e}^{\text{i}2lk_{c}x}+\text{c.c.}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{X}_{2n+1}=\mathop{\sum }_{l=0}^{n}\boldsymbol{Z}_{2n+1,2l+1}(z)\text{e}^{\text{i}(2l+1)k_{c}x}+\text{c.c.}, & \displaystyle\end{eqnarray}$$
the vertical normal mode
$\boldsymbol{Z}_{n,l}=(P_{n,l}(z);U_{n,l}(z);W_{n,l}(z);\unicode[STIX]{x1D6E9}_{n,l}(z))^{\text{T}}$
being indexed with the order
$n$
of the solution and harmonic number
$l$
in the
$x$
dependence. It also appears recursively that
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{2n}=-{\displaystyle \frac{\displaystyle \mathop{\sum }_{l=1}^{2n}\langle \boldsymbol{X}_{1}|\boldsymbol{N}(\boldsymbol{X}_{l},\boldsymbol{X}_{2n+1-l})\rangle +\displaystyle \mathop{\sum }_{l=1}^{n-1}\mathit{Ra}_{2l}\langle \boldsymbol{X}_{1}|\unicode[STIX]{x1D648}\boldsymbol{X}_{2(n-l)+1}\rangle }{\langle \boldsymbol{X}_{1}|\unicode[STIX]{x1D648}\boldsymbol{X}_{1}\rangle }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{2n+1}=0. & \displaystyle\end{eqnarray}$$
This is true for orders 1 and 2, as explained above and, assuming it holds up to degrees
$2n-1$
and
$2n$
, the expressions for degrees
$2n+1$
and
$2n+2$
can be predicted from (4.30). First, equation (4.30) of order
$2n+1$
includes on the right-hand side only terms up to degree
$2n$
and can be used to predict the form of
$\boldsymbol{X}_{2n+1}$
. Each term of the form
$\boldsymbol{N}(\boldsymbol{X}_{l},\boldsymbol{X}_{2n+1-l})$
contains only odd powers of
$\text{e}^{\text{i}k_{c}x}$
since it is composed of products of even (respectively odd) and odd (respectively even) polynomials of
$\text{e}^{\text{i}k_{c}x}$
for
$l$
even (respectively odd). Each term of the form
$\mathit{Ra}_{l}\unicode[STIX]{x1D648}\boldsymbol{X}_{2n+1-l}$
is either null for
$l$
odd or an odd polynomial of
$\text{e}^{\text{i}k_{c}x}$
for
$l$
even. Summing up, the right-hand side of the equation being an odd polynomial of
$\text{e}^{\text{i}k_{c}x}$
, the solution to the equation is of the form (4.36).
Taking the dot product of (4.30) of order
$2n+2$
with
$\boldsymbol{X}_{1}$
and using the Hermitian character of
$\unicode[STIX]{x1D647}_{c}$
provides the equation for
$\mathit{Ra}_{2n+1}$
. Starting first with the last term on the right-hand side, all the terms in the sum except the one in
$\mathit{Ra}_{2n+1}$
drop out either because
$\mathit{Ra}_{l}$
is null for
$l$
odd or because the dot product
$\langle \boldsymbol{X}_{1}|\unicode[STIX]{x1D648}\boldsymbol{X}_{2n+2-l}\rangle =0$
for
$l$
even, since
$\boldsymbol{X}_{2n+2-l}$
then contains only even powers of
$\text{e}^{\text{i}k_{c}x}$
. We are left with
$\mathit{Ra}_{2n+1}\langle \boldsymbol{X}_{1}|\unicode[STIX]{x1D648}\boldsymbol{X}_{1}\rangle$
. Considering the first sum on the right-hand side, each term
$\boldsymbol{N}(\boldsymbol{X}_{l},\boldsymbol{X}_{2n+2-l})$
is an even polynomial of
$\text{e}^{\text{i}k_{c}x}$
, as the product of either two even polynomials (for
$l$
even) or two odd polynomials (for
$l$
odd). Therefore, each of these terms is orthogonal to
$\boldsymbol{X}_{1}$
and
$\mathit{Ra}_{2n+1}=0$
. The same equation (4.30)
$_{2n+2}$
contains only even powers of
$\text{e}^{\text{i}k_{c}x}$
on the right-hand side, and this justifies equation (4.35) for the order
$2n+2$
.
Finally, (4.37)
$_{2n+2}$
is obtained by simply taking the dot product of (4.30)
$_{2n+3}$
with
$\boldsymbol{X}_{1}$
.
An important diagnostic for convection is the heat transfer efficiency measured by the dimensionless mean heat flux density, the Nusselt number
$\mathit{Nu}$
. Since the temperature is uniform on each horizontal boundary and the average vertical velocity is null for the deforming mode considered here, the advective heat transfer across the horizontal boundaries is null. Therefore, the Nusselt number can easily be computed by taking the vertical derivative of the temperature at either boundary. In the Fourier decomposition used for the nonlinear analysis, only the zeroth order terms in
$\text{e}^{\text{i}k_{c}x}$
contribute to the horizontal average, and they only appear in terms that are even in the
$\unicode[STIX]{x1D716}$
development. Restricting ourselves here to an order-2 development, the Nusselt number can be computed as
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Nu}=1-\unicode[STIX]{x1D716}^{2}D\unicode[STIX]{x1D703}_{2,0}\left({\displaystyle \frac{1}{2}}\right)=1-D\unicode[STIX]{x1D703}_{2,0}\left({\displaystyle \frac{1}{2}}\right){\displaystyle \frac{\mathit{Ra}_{c}}{\mathit{Ra}_{2}}}{\displaystyle \frac{\mathit{Ra}-\mathit{Ra}_{c}}{\mathit{Ra}_{c}}}, & \displaystyle\end{eqnarray}$$
where (4.34) was used. This equation shows the classical result that the convective heat flow,
$\mathit{Nu}-1$
, increases linearly with the reduced Rayleigh number
$\unicode[STIX]{x1D700}=(\mathit{Ra}-\mathit{Ra}_{c})/\mathit{Ra}_{c}$
for small values of
$\unicode[STIX]{x1D700}$
and the determination of the coefficient of proportionality,
$A$
, is the main goal of the weakly nonlinear analysis presented here. Note that
$\boldsymbol{N}(\boldsymbol{X}_{2},\boldsymbol{X}_{1})$
and
$\boldsymbol{N}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})$
have a non-zero component only along the
$\unicode[STIX]{x1D703}$
space (4.23) so that, because of our definition of the dot product (appendix A) and using (4.33),
$\mathit{Ra}_{2}$
is proportional to
$\mathit{Ra}_{c}$
.
The procedure just outlined can be applied to the case with classical boundary conditions. In particular, for free-slip non-penetrating boundary conditions, the problem can be solved analytically (Malkus & Veronis Reference Malkus and Veronis1958; Manneville Reference Manneville2004). Starting with the vertical velocity in the critical mode as
$w_{1}=\sin kx\cos \unicode[STIX]{x03C0}z$
, one gets
$\unicode[STIX]{x1D703}_{1}=(\unicode[STIX]{x03C0}^{2}+k^{2})^{-1}\sin kx\cos \unicode[STIX]{x03C0}z$
,
$\mathit{Ra}_{c}=(\unicode[STIX]{x03C0}^{2}+k^{2})^{3}/k^{2}$
,
$w_{2}=0$
,
$\unicode[STIX]{x1D703}_{2}=(8\unicode[STIX]{x03C0}(\unicode[STIX]{x03C0}^{2}+k^{2}))^{-1}\sin 2\unicode[STIX]{x03C0}z$
and
$\mathit{Ra}_{2}=(\unicode[STIX]{x03C0}^{2}+k^{2})^{2}/8k^{2}$
. This gives
$A=-D\unicode[STIX]{x1D703}_{2,0}(1/2)\mathit{Ra}_{c}/\mathit{Ra}_{2}=2$
.
Heat flux coefficient as a function of the phase-change numbers, equal to each other (a), and Nusselt number as a function of Rayleigh number for different values of
$\unicode[STIX]{x1D6F7}^{\pm }$
(b). In (a), the solid line gives the limit of two non-penetrating boundaries while the dashed line represents the first-order development obtained for
$\unicode[STIX]{x1D6F7}\rightarrow 0$
(4.44).
Similarly, the low-
$\unicode[STIX]{x1D6F7}$
expansion of the linear mode, equations (4.14)–(4.19), can be used to compute the behaviour of coefficient
$A$
at low
$\unicode[STIX]{x1D6F7}$
values. We choose
$\unicode[STIX]{x1D6E9}_{max}=1/16$
to have a normalisation consistent with the one above (the amplitude of
$\boldsymbol{X}_{1}$
is not defined by the linear problem and changing its normalisation, say by multiplying it by a factor
$a$
, leads to
$\boldsymbol{X}_{2}$
and
$\mathit{Ra}_{2}$
multiplied by
$a^{2}$
, so that by virtue of (4.34), the total solution
$\boldsymbol{X}$
is unchanged) and the solution at order 2 is searched for in the form of
$z$
polynomials, and we get, to order 1 in
$\unicode[STIX]{x1D6F7}$
,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{2}=-{\displaystyle \frac{z}{48}}\left(z^{2}-{\displaystyle \frac{1}{4}}\right)\left[1+\left(1-{\displaystyle \frac{\unicode[STIX]{x1D6F7}}{64}}\right)\cos 2k_{c}x\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{2}=-{\displaystyle \frac{\sqrt{\unicode[STIX]{x1D6F7}}}{192\sqrt{2}}}\sin 2k_{c}x, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle w_{2}={\displaystyle \frac{z\unicode[STIX]{x1D6F7}}{256}}\cos 2k_{c}x, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}_{2}={\displaystyle \frac{1}{320}}-{\displaystyle \frac{43\unicode[STIX]{x1D6F7}}{430\,080}}. & \displaystyle\end{eqnarray}$$
The heat flux coefficient is then, to order 1 in
$\unicode[STIX]{x1D6F7}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle A={\displaystyle \frac{4480}{1344-43\unicode[STIX]{x1D6F7}}}. & \displaystyle\end{eqnarray}$$
Figure 11(a) represents the value of the heat flux coefficient
$A$
as a function of
$\unicode[STIX]{x1D6F7}$
obtained using the Chebyshev-collocation approach described above (solid circles, see appendix B for details on the calculation of nonlinear terms) and the two limiting cases of
$\unicode[STIX]{x1D6F7}\rightarrow \infty$
(solid line) and
$\unicode[STIX]{x1D6F7}\rightarrow 0$
(dashed line), which shows a good match.
The heat flux coefficient
$A$
, which equals 2 for classical non-penetrating boundaries, tends to
$10/3$
when
$\unicode[STIX]{x1D6F7}\rightarrow 0$
. This
${\sim}50\,\%$
increase makes the Nusselt number increase when
$\unicode[STIX]{x1D6F7}$
tends to zero, but the main effect comes from the decrease of the critical Rayleigh number as
${\sim}24\unicode[STIX]{x1D6F7}$
, which makes the slope
$\text{d}\mathit{Nu}/\text{d}\mathit{Ra}$
go to infinity as
${\sim}5/36\unicode[STIX]{x1D6F7}$
. This is illustrated on figure 11(b), which shows the
$\mathit{Nu}-\mathit{Ra}$
relationship derived from this analysis for different values of
$\unicode[STIX]{x1D6F7}$
. The heat transfer efficiency is greatly increased by decreasing
$\unicode[STIX]{x1D6F7}$
for two reasons. First, it makes the critical Rayleigh number decrease so that convection starts with a lower Rayleigh number. Second, the rate at which the Nusselt number increases with
$\mathit{Ra}$
above its critical value is also drastically increased when
$\unicode[STIX]{x1D6F7}$
is decreased.
First unstable mode when only the bottom boundary is a phase-change interface, with
$\unicode[STIX]{x1D6F7}^{-}=10$
(a) and
$\unicode[STIX]{x1D6F7}^{-}=10^{-2}$
(b). The temperature anomaly compared to the conduction solution is represented in colours and streamlines have a thickness proportional to the relative norm of the velocity.
5 Solutions with only one phase-change boundary
Let us now consider the case when only one boundary is a liquid–solid phase change, the other one being subject to a non-penetrating condition. With the plane layer geometry considered here, the situation with the upper boundary having a phase change is symmetrical to the one with a lower boundary having a phase change. The latter is considered here since it applies to the dynamics of the icy shells of some satellites of giant planets (Čadek et al. Reference Čadek, Tobie, Van Hoolst, Massé, Choblet, Lefèvre, Mitri, Baland, Běhounková, Bourgeois and Trinh2016) and possibly to Earth’s mantle for a large part of its history (Labrosse et al. Reference Labrosse, Hernlund and Coltice2007).
Critical Rayleigh number (a) and wavenumber (b) as function of the phase-change number for the bottom boundary
$\unicode[STIX]{x1D6F7}^{-}$
, the top boundary having a non-penetrating condition. The dash-dotted lines represent the classical values obtained for two non-penetrating conditions, for reference.
The analysis is done in the same way as for the case with a phase change at both boundaries. Figure 12 shows examples of the first unstable mode for two different values of
$\unicode[STIX]{x1D6F7}^{-}$
. The upper one shows that when
$\unicode[STIX]{x1D6F7}^{-}=10$
, the convection geometry is not very different from that with a non-penetrating condition (hereafter ‘the classical situation’) but the streamlines are slightly open at the bottom. The horizontal wavelength at onset,
$\unicode[STIX]{x1D706}_{c}=3.57$
, is larger than the one for the classical situation (
$\unicode[STIX]{x1D706}_{c}=2\sqrt{2}$
) and the critical Rayleigh number is smaller (
$\mathit{Ra}_{c}=352$
). For
$\unicode[STIX]{x1D6F7}^{-}=10^{-2}$
, the streamlines are almost normal to the bottom boundary and the wavelength
$\unicode[STIX]{x1D706}_{c}=5$
is approximately twice the classical one, as if the solution was the upper half of a classical convective domain. However, the boundary condition imposed for temperature at the bottom is different from what would be obtained in that case, and the critical Rayleigh number,
$\mathit{Ra}_{c}=153$
, is approximately a quarter of the classical one. This can be understood in a heuristic way: the Rayleigh number can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Ra}={\displaystyle \frac{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}}{\unicode[STIX]{x1D70F}_{c}^{2}}}={\displaystyle \frac{\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}Tg}{d}}{\displaystyle \frac{d^{2}}{\unicode[STIX]{x1D708}}}{\displaystyle \frac{d^{2}}{\unicode[STIX]{x1D705}}}, & \displaystyle\end{eqnarray}$$
with
$\unicode[STIX]{x1D70F}_{c}$
the convective time scale associated with acceleration due to gravity,
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}$
the viscous time scale and
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}$
the thermal diffusion time scale. Compared to the classical situation, we have the same imposed temperature gradient, hence the same
$\unicode[STIX]{x1D70F}_{c}$
. Similarly, diffusion happens on the same vertical length scale and we have the same
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}$
. On the other hand, the bottom boundary imposes no limit to vertical flow and the viscous deformation is distributed over a vertical distance twice the thickness of the layer, which increases the effective viscous time scale by a factor of four. Therefore, the Rayleigh number imposed here is equivalent to a value four times larger than in the classical situation.
Figure 13 shows the variation of the critical Rayleigh number (panel a) and wavenumber (panel b) as a function of
$\unicode[STIX]{x1D6F7}^{-}$
, and one can see that both tend to a finite value when
$\unicode[STIX]{x1D6F7}^{-}\rightarrow 0$
. The mode obtained for
$\unicode[STIX]{x1D6F7}^{-}=10^{-2}$
is close to that limit. Contrary to the situation with a phase change at both boundaries, the presence of a non-penetrating boundary condition implies that some deformation is always needed for convection to occur, which makes viscosity still a limiting factor at vanishing values of
$\unicode[STIX]{x1D6F7}^{-}$
.
Heat flux coefficient as a function of the bottom phase-change number
$\unicode[STIX]{x1D6F7}^{-}$
, the top boundary being non-penetrative (a), and Nusselt number as a function of Rayleigh number for different values of
$\unicode[STIX]{x1D6F7}^{-}$
(b).
Considering now the weakly nonlinear analysis results, figure 14(a) shows that the heat flux coefficient for only one phase-change boundary condition tends to slightly greater than 1; that is, approximately half that for the case for both non-penetrative boundaries. Combining that with a critical Rayleigh number that is approximately four times smaller makes
$\text{d}\mathit{Nu}/\text{d}\mathit{Ra}$
approximately twice that for the classical situation. Therefore, the efficiency of heat transfer is improved compared to the classical case, both because convection starts for a smaller Rayleigh number and because the rate of variation of the Nusselt number with
$\mathit{Ra}$
is approximately twice as large. This is illustrated on figure 14(b).
Mean temperature coefficient (
$B$
defined in (5.2)) as function of the bottom phase-change parameter
$\unicode[STIX]{x1D6F7}^{-}$
(a) and mean temperature as function of
$\mathit{Ra}$
for different values of
$\unicode[STIX]{x1D6F7}^{-}$
(b). The range of
$\mathit{Ra}$
values explored is the same as that used for figure 14.
In contrast to the case with both boundaries being a phase change with equal values of
$\unicode[STIX]{x1D6F7}$
, the case discussed in this section breaks the symmetry around the
$z=0$
plane. In particular, this means that the mean temperature in the domain is not equal to the average of both boundaries,
$\langle T\rangle \neq 1/2$
in dimensionless form. As for the Nusselt number (4.39), a contribution from all even orders in
$\unicode[STIX]{x1D716}$
is expected, and to the leading order explored here,
$$\begin{eqnarray}\displaystyle & \displaystyle \langle T\rangle ={\displaystyle \frac{1}{2}}+\langle \unicode[STIX]{x1D703}_{2,0}\rangle {\displaystyle \frac{\mathit{Ra}_{c}}{\mathit{Ra}_{2}}}{\displaystyle \frac{\mathit{Ra}-\mathit{Ra}_{c}}{\mathit{Ra}_{c}}}\equiv {\displaystyle \frac{1}{2}}+B{\displaystyle \frac{\mathit{Ra}-\mathit{Ra}_{c}}{\mathit{Ra}_{c}}}. & \displaystyle\end{eqnarray}$$
The coefficient
$B$
defined above is computed exactly for the case of both non-penetrating boundaries, and as expected found to be null. Figure 15(a) shows the evolution of this coefficient as a function of
$\unicode[STIX]{x1D6F7}^{-}$
. One can see that it tends to a finite positive value in the limit
$\unicode[STIX]{x1D6F7}^{-}\rightarrow 0$
. Therefore, for small values of
$\unicode[STIX]{x1D6F7}^{-}$
, the average temperature is expected to be larger than
$1/2$
(figure 15
b). For the same range of Rayleigh number as explored in figure 14, figure 15(b) also shows the evolution of the mean temperature at the leading order given by (5.2). For low values of
$\unicode[STIX]{x1D6F7}^{-}$
, the mean temperature increases rapidly with Rayleigh number.
Finite-amplitude solution for
$\unicode[STIX]{x1D6F7}^{-}=10^{-2}$
,
$\unicode[STIX]{x1D716}=5.58$
and a non-penetrating boundary condition at the top.
The asymmetry of the mean temperature for low values of
$\unicode[STIX]{x1D6F7}^{-}$
is also expressed in the finite-amplitude solution that can be plotted for a given value of
$\unicode[STIX]{x1D716}$
. The range of validity of such solutions as a function of
$\unicode[STIX]{x1D716}$
depends on the order of the development. Computing the solution only up to order 3 in
$\unicode[STIX]{x1D716}$
, we restrict ourselves to small values of this number and figure 16 shows the result for
$\unicode[STIX]{x1D716}=5.58$
corresponding to
$\mathit{Nu}=1.5$
. This shows that the down-welling current is more focused than the up-welling one. This situation is similar to the case of volumetrically heated convection (e.g. Parmentier & Sotin Reference Parmentier and Sotin2000), which is not the case here. Preliminary direct numerical simulations confirm this behaviour, but the full exploration of this question goes beyond the scope of the present paper.
6 Conclusion
In the context of the dynamics of planetary mantles, convection can happen in solid shells adjacent to liquid layers. The viscous stress in the solid builds up a topography of the interface between the solid and liquid layers. In the absence of mechanisms to erase topography, its buoyancy equilibrates the viscous stress, which effectively enforces a non-penetrating boundary condition. On the other hand, if the topography can be suppressed by melting and freezing at the interface at a faster pace than its building process, the vertical velocity is not required to be null at the interface. The non-penetrating boundary condition is then replaced by a relationship between the normal velocity, its normal gradient and pressure (2.16), and involving a dimensionless phase-change number,
$\unicode[STIX]{x1D6F7}$
, the ratio of the phase-change time scale to the viscous time scale (2.15). When this number is large, we recover the classical non-penetrating condition, while the limit of low
$\unicode[STIX]{x1D6F7}$
allows a large flow through the boundary.
When both boundaries are characterised by a low
$\unicode[STIX]{x1D6F7}$
number, a translating, non-deforming, mode of convection is possible and competes with a deforming mode with wavenumber
$k$
that decreases as
$\sqrt{\unicode[STIX]{x1D6F7}}$
, and therefore ressembles translation with alternating up- and downward directions. The critical Rayleigh number for the onset of the deforming mode is slightly below that of the translation mode,
$\mathit{Ra}=24\unicode[STIX]{x1D6F7}$
, but the latter is found to be stable against a deforming instability when the Rayleigh number is
${\sim}\unicode[STIX]{x1D6F7}^{2}$
above the critical value. It is therefore likely to dominate when both boundaries are characterised by low values of
$\unicode[STIX]{x1D6F7}$
. In both translating and deforming modes of convection, the heat transfer efficiency, the Nusselt number, is found to increase strongly with Rayleigh number at small values of
$\unicode[STIX]{x1D6F7}$
.
When only one boundary is a phase-change interface with a low value of
$\unicode[STIX]{x1D6F7}$
, the wavenumber is approximately half and the critical Rayleigh number is approximately a quarter those of the corresponding values for the classical non-penetrating boundary condition. Close to onset, a weakly nonlinear analysis shows that the Nusselt number varies linearly with the Rayleigh number with a slope that is approximately twice that for both non-penetrating boundary conditions. The average temperature is also found to increase strongly with Rayleigh number and the flow geometry is strongly affected, with down-welling currents more focused than up-welling ones.
Overall, having the possibility of melting and freezing across one or both horizontal boundaries of an infinite Prandtl number fluid makes convection much easier (i.e. the critical Rayleigh number is strongly reduced), the preferred horizontal wavelength much larger and heat transfer much stronger, with important potential implications for planetary dynamics.
Acknowledgements
We are grateful to three anonymous reviewers and editor G. Worster for comments that pushed us to significantly clarify our paper. This research has been funded by the French Agence Nationale de la Recherche under the grant number ANR-15-CE31-0018-01, MaCoMaOc.
Appendix A. Self-adjointness of operator
$\unicode[STIX]{x1D647}_{c}$
Using a Fourier decomposition for the horizontal decomposition,
$\unicode[STIX]{x1D647}_{c}$
simply reads as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D647}_{c}=\left(\begin{array}{@{}cccc@{}}0 & \text{i}k & \text{D} & 0\\ -\text{i}k & \text{D}^{2}-k^{2} & 0 & 0\\ -\text{D} & 0 & \text{D}^{2}-k^{2} & \mathit{Ra}_{c}\\ 0 & 0 & 1 & \text{D}^{2}-k^{2},\end{array}\right) & \displaystyle\end{eqnarray}$$
where the time derivative has been omitted since the linear instability is found to be stationary. In a linear stability analysis, adding a growth rate
$\unicode[STIX]{x1D70E}$
on the diagonal of the matrix would not alter the adjoint calculation, as will appear below. The boundary conditions are given by (4.5) to (4.7). In the calculation of the dot product, the
$\unicode[STIX]{x1D703}$
part is given
$\mathit{Ra}_{c}$
as weight and the horizontal integral can be factored out:
$$\begin{eqnarray}\displaystyle \displaystyle \langle \boldsymbol{X}_{2}|\unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{1}\rangle & = & \displaystyle \int \text{e}^{\text{i}(k_{2}-k_{1})}\text{d}x\left[\int _{-1/2}^{1/2}\bar{P}_{2}(\text{i}kU_{1}+\text{D}W_{1})\,\text{d}z\right.\nonumber\\ \displaystyle & & \displaystyle +\,\int _{-1/2}^{1/2}\bar{U}_{2}(-\text{i}kP_{1}+(\text{D}^{2}-k^{2})U_{1})\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle +\,\int _{-1/2}^{1/2}\bar{W}_{2}(-\text{D}P_{1}+(\text{D}^{2}-k^{2})W_{1}+\mathit{Ra}_{c}\unicode[STIX]{x1D6E9}_{1})\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle \left.+\,\mathit{Ra}_{c}\int _{-1/2}^{1/2}\bar{\unicode[STIX]{x1D6E9}}_{2}(W_{1}+(\text{D}^{2}-k^{2})\unicode[STIX]{x1D6E9}_{1})\,\text{d}z\right],\end{eqnarray}$$
where the overbar means complex conjugate. Since the
$x$
part poses no difficulty, we only consider the
$z$
part, which we denote as
$\langle \bullet |\bullet \rangle _{z}$
. Reordering the different integrals in (A 2) so that terms of
$\boldsymbol{X}_{1}$
are factored out and performing integration by parts on each term including
$\text{D}$
, we get
$$\begin{eqnarray}\displaystyle \displaystyle \langle \boldsymbol{X}_{2}|\unicode[STIX]{x1D647}_{c}\boldsymbol{X}_{1}\rangle _{z} & = & \displaystyle \int _{-1/2}^{1/2}(-\text{i}k\bar{U}_{2}+\text{D}\bar{W}_{2})P_{1}\,\text{d}z+\int _{-1/2}^{1/2}(\text{i}k\bar{P}_{2}+(\text{D}^{2}-k^{2})\bar{U}_{2})U_{1}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle +\,\int _{-1/2}^{1/2}(\text{D}\bar{P}_{2}+(\text{D}^{2}-k^{2})\bar{W}_{1}+\mathit{Ra}_{c}\bar{\unicode[STIX]{x1D6E9}}_{2})W_{1}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle +\,\mathit{Ra}_{c}\int _{-1/2}^{1/2}(\bar{W}_{2}+(\text{D}^{2}-k^{2})\bar{\unicode[STIX]{x1D6E9}}_{2})\unicode[STIX]{x1D6E9}_{1}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle +\,[\bar{P}_{2}W_{1}]_{-1/2}^{1/2}+[\bar{U}_{2}\text{D}U_{1}]_{-1/2}^{1/2}-[U_{1}\text{D}\bar{U}_{2}]_{-1/2}^{1/2}-[\bar{W}_{2}P_{1}]_{-1/2}^{1/2}\nonumber\\ \displaystyle & & \displaystyle +\,[\bar{W}_{2}\text{D}W_{1}]_{-1/2}^{1/2}-[W_{1}\text{D}\bar{W}_{2}]_{-1/2}^{1/2}+\mathit{Ra}([\bar{\unicode[STIX]{x1D6E9}}_{2}\text{D}\unicode[STIX]{x1D6E9}_{1}]_{-1/2}^{1/2}-[\unicode[STIX]{x1D6E9}_{1}\text{D}\bar{\unicode[STIX]{x1D6E9}}_{2}]_{-1/2}^{1/2}).\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The integral part shows that the adjoint linear system is the same as the direct one, with
$\unicode[STIX]{x1D647}_{c}$
as operator. The boundary conditions are such as to allow suppression of all the boundary values in (A 3). The boundary conditions (4.5) to (4.7) are applied to
$\boldsymbol{X}_{1}$
to remove
$\unicode[STIX]{x1D6E9}_{1}(\pm 1/2)$
and replace
$\text{D}U_{1}$
and
$P_{1}$
. In addition, the mass conservation equation applied to
$\boldsymbol{X}_{2}$
allows one to replace
$\text{D}W_{2}$
. Factorising
$W_{1}$
,
$U_{1}$
and
$\unicode[STIX]{x1D6E9}_{1}$
gives for the boundary conditions
$$\begin{eqnarray}\displaystyle & \displaystyle [W_{1}(-\bar{P}_{2}\pm \unicode[STIX]{x1D6F7}^{\pm }\bar{W}_{2}+2\text{D}\bar{W}_{2})]_{-1/2}^{1/2}+[U_{1}(-\text{i}k\bar{W}_{2}+\text{D}\bar{U}_{2})]_{-1/2}^{1/2}-\mathit{Ra}[\bar{\unicode[STIX]{x1D6E9}}_{2}\text{D}\unicode[STIX]{x1D6E9}_{1}]_{-1/2}^{1/2}=0. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Since
$W_{1}$
,
$U_{1}$
and
$\text{D}\unicode[STIX]{x1D6E9}_{1}$
can take arbitrary values on the boundaries, the differences can only be eliminated in a general manner by setting all their coefficients to 0, which gives the boundary conditions for the adjoint:
$$\begin{eqnarray}\displaystyle & \displaystyle \text{D}U_{2}+\text{i}kW_{2}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \pm \unicode[STIX]{x1D6F7}^{\pm }W_{2}+2\text{D}W_{2}-P_{2}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}_{2}=0. & \displaystyle\end{eqnarray}$$
The adjoint problem is therefore identical to the direct one. Among other implications, all eigenvalues of
$L_{c}$
must be real, which is consistent with our numerical findings.
Appendix B. Expression of the nonlinear terms
Computation of the nonlinear term
$\boldsymbol{N}(\boldsymbol{X}_{n},\boldsymbol{X}_{m})$
(4.23) is the trickiest part of the procedure explained in § 4.2 and deserves some details provided here. First of all, it contains only a
$\unicode[STIX]{x1D723}$
component, referred to as
$\boldsymbol{N}(\boldsymbol{X}_{n},\boldsymbol{X}_{m})_{\unicode[STIX]{x1D6E9}}$
. To compute it, one needs first to decompose indices
$n$
and
$m$
as
$$\begin{eqnarray}\displaystyle & \displaystyle n=2p+q\quad \text{with }p=\left\lfloor {\displaystyle \frac{n}{2}}\right\rfloor , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle m=2r+s\quad \text{with }r=\left\lfloor {\displaystyle \frac{m}{2}}\right\rfloor , & \displaystyle\end{eqnarray}$$
where
$\lfloor \,\rfloor$
denotes the floor function. In computing
$\boldsymbol{N}(\boldsymbol{X}_{n},\boldsymbol{X}_{m})_{\unicode[STIX]{x1D6E9}}$
, one needs to account for the full (i.e. real) expression of
$\boldsymbol{X}_{n}$
and
$\boldsymbol{X}_{m}$
including the complex conjugate. They write
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{X}_{n}=\mathop{\sum }_{l_{1}=0}^{p}\boldsymbol{Z}_{n,2l_{1}+q}(z)\text{e}^{\text{i}(2l_{1}+q)kx}+\text{c.c.}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{X}_{m}=\mathop{\sum }_{l_{2}=0}^{r}\boldsymbol{Z}_{m,2l_{1}+s}(z)\text{e}^{\text{i}(2l_{2}+s)kx}+\text{c.c.} & \displaystyle\end{eqnarray}$$
Using (4.23), we get
$$\begin{eqnarray}\displaystyle \boldsymbol{N}(\boldsymbol{X}_{n},\boldsymbol{X}_{m})_{\unicode[STIX]{x1D6E9}} & = & \displaystyle \mathop{\sum }_{l_{1}=0}^{p}\mathop{\sum }_{l_{2}=0}^{q}\big\{\!\big[\text{i}(2l_{2}+s)kU_{n,2l_{1}+q}\unicode[STIX]{x1D6E9}_{m,2l_{2}+s}\nonumber\\ \displaystyle & & \displaystyle +\,W_{n,2l_{1}+q}D\unicode[STIX]{x1D6E9}_{m,2l_{2}+s}\big]\text{e}^{\text{i}[2(l_{1}+l_{2})+q+s]kx}\nonumber\\ \displaystyle & & \displaystyle +\,\big[-\text{i}(2l_{2}+s)kU_{n,2l_{1}+q}\bar{\unicode[STIX]{x1D6E9}}_{m,2l_{2}+s}\nonumber\\ \displaystyle & & \displaystyle +\,W_{n,2l_{1}+q}D\bar{\unicode[STIX]{x1D6E9}}_{m,2l_{2}+s}\big]\text{e}^{\text{i}(2(l_{1}-l_{2})+q-s)kx}\big\}+\text{c.c.}\end{eqnarray}$$
The harmonics of the first term are always positive while those of the second term can be negative. Either way, each term has its complex conjugate and we solve only for the positive or null harmonics, the rest of the solution simply being obtained as the conjugate of the computed part.
