3.1 Phase-field approach for the interface

In order to solve the previous dimensionless equations without having to impose the internal boundary conditions related to the interface, we introduce the continuous phase-field or order parameter $\unicode[STIX]{x1D719}(x,z,t)$ such that $\unicode[STIX]{x1D719}=0$ in the solid phase and $\unicode[STIX]{x1D719}=1$ in the liquid. A thin interface of finite width in which $\unicode[STIX]{x1D719}$ takes values between zero and unity exists between the pure solid and liquid phases. Writing an evolution equation for the phase-field parameter can be done in several ways. The simpler derivation, which we briefly explain here, is the geometrical approach described in Beckermann et al. (Reference Beckermann, Diepers, Steinbach, Karma and Tong1999) and starts from the Gibbs–Thomson effect

(3.1) $$\begin{eqnarray}\displaystyle \frac{v_{n}}{\unicode[STIX]{x1D707}}=T_{M}-T-\frac{\unicode[STIX]{x1D70E}_{s}T_{M}}{{\mathcal{L}}}\unicode[STIX]{x1D705}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is the mobility, $\unicode[STIX]{x1D70E}_{s}$ the surface tension, $\unicode[STIX]{x1D705}$ the mean curvature of the front and $v_{n}$ the normal velocity of the interface between the solid and the liquid phases. Although $\unicode[STIX]{x1D719}$ represents a finite-thickness interface, the normal velocity of the front can be related to the time evolution of $\unicode[STIX]{x1D719}$ at a fixed value (for instance $\unicode[STIX]{x1D719}=1/2$ ), through the equation

(3.2) $$\begin{eqnarray}\displaystyle v_{n}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}t}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}. & & \displaystyle\end{eqnarray}$$

Moreover, the curvature of the front can be computed in terms of $\unicode[STIX]{x1D719}$ through

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|}\right)_{\unicode[STIX]{x1D719}=1/2}. & & \displaystyle\end{eqnarray}$$

Substituting (3.2) and (3.3) into (3.1), we obtain an evolution equation for $\unicode[STIX]{x1D719}$ , in which the right-hand side depends only on $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}$ . However, this equation does not have a unique stationary solution. Therefore, the profile for $\unicode[STIX]{x1D719}$ has to be specified, and this point is motivated by thermodynamics considerations.

This leads us to the second approach for deriving the evolution equation for $\unicode[STIX]{x1D719}$ , based on thermodynamics and described in detail in Wang et al. (Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993) among others (Penrose & Fife Reference Penrose and Fife1990; Karma & Rappel Reference Karma and Rappel1996). The entropy of a given volume $V$ is represented by the functional

(3.4) $$\begin{eqnarray}\displaystyle {\mathcal{S}}=\int _{V}\left[s-\frac{\unicode[STIX]{x1D6FF}^{2}}{2}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})^{2}\right]\,\text{d}V, & & \displaystyle\end{eqnarray}$$

where $s(e,\unicode[STIX]{x1D719})$ is the entropy density, $e$ is the internal energy density, $\unicode[STIX]{x1D719}$ the phase field and $\unicode[STIX]{x1D6FF}$ a constant. The second term in the right-hand side of equation (3.4) is analogous to the Landau–Ginzburg gradient term in the free energy and is accounting for contributions from the liquid–solid interface. In order to ensure that the local entropy production is positive (Wang et al. Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993), the phase field must evolve according to

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}=\left.\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\right|_{e}+\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}$ is a positive constant. Following the thermodynamically consistent derivation of Wang et al. (Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993), this leads to the following dimensional phase-field equation:

(3.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}+Q(T)\frac{\text{d}p(\unicode[STIX]{x1D719})}{\text{d}\unicode[STIX]{x1D719}}-\frac{1}{4a}\frac{\text{d}g(\unicode[STIX]{x1D719})}{\text{d}\unicode[STIX]{x1D719}}, & & \displaystyle\end{eqnarray}$$

where $Q(T)$ is defined as

(3.7) $$\begin{eqnarray}\displaystyle Q(T)=\int _{T_{M}}^{T}\frac{{\mathcal{L}}(\unicode[STIX]{x1D701})}{\unicode[STIX]{x1D701}^{2}}\,\text{d}\unicode[STIX]{x1D701}. & & \displaystyle\end{eqnarray}$$

In the following, we assume that the latent heat ${\mathcal{L}}$ does not depend on temperature and that the temperature close to the interface is always approximately the melting temperature $T_{M}$ , i.e. $|T-T_{M}|\ll T_{M}$ , so that (3.7) can be simplified to

(3.8) $$\begin{eqnarray}\displaystyle Q(T)\approx \frac{{\mathcal{L}}}{T_{M}^{2}}(T-T_{M}). & & \displaystyle\end{eqnarray}$$

Note that the validity of this simplification can be questionable in our case since thermal boundary layers will develop close to the interface. We nevertheless checked its impact on our results by comparing the original function defined by (3.7) to its simplified version (3.8), and found no significant differences.

The two functions $p(\unicode[STIX]{x1D719})$ and $g(\unicode[STIX]{x1D719})$ must be prescribed in order to close the model. While several choices exist in the literature, we use the prescription of Wang et al. (Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993) which ensures that the solid and liquid phases correspond to $\unicode[STIX]{x1D719}=0$ and $\unicode[STIX]{x1D719}=1$ , irrespective of the temperature distribution across both phases:

(3.9) $$\begin{eqnarray}\displaystyle g(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D719}^{2}(1-\unicode[STIX]{x1D719})^{2} & & \displaystyle\end{eqnarray}$$

and

(3.10) $$\begin{eqnarray}\displaystyle p(\unicode[STIX]{x1D719})=\frac{\displaystyle \int _{0}^{\unicode[STIX]{x1D719}}g(\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D709}}{\displaystyle \int _{0}^{1}g(\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D709}}=\unicode[STIX]{x1D719}^{3}(10-15\unicode[STIX]{x1D719}+6\unicode[STIX]{x1D719}^{2}). & & \displaystyle\end{eqnarray}$$

The function $g(\unicode[STIX]{x1D719})$ corresponds to a double well and ensures that the phase field is either equal to $0$ or $1$ everywhere except close to the liquid–solid interface where the phase change occurs. The positive constant $a$ in (3.6) is related to the amplitude of the potential barrier between the two equilibria. The function $p(\unicode[STIX]{x1D719})$ ensures a continuous transition between each extremum value of $\unicode[STIX]{x1D719}$ . Note that in a steady one-dimensional configuration, and assuming that $T=T_{M}$ , equation (3.6) leads to a simple analytical profile for the phase variable around the interface located at $x=x_{i}$ given by

(3.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}(x)=\frac{1}{2}\left[1-\tanh \left(\frac{x-x_{i}}{2\sqrt{2a}\unicode[STIX]{x1D6FF}}\right)\right], & & \displaystyle\end{eqnarray}$$

assuming that $\unicode[STIX]{x1D719}=1$ as $x\rightarrow -\infty$ and $\unicode[STIX]{x1D719}=0$ as $x\rightarrow +\infty$ . The diffuse interface therefore has a characteristic thickness equal to $\unicode[STIX]{x1D6FF}\sqrt{a}$ .

The corresponding dimensional temperature equation is given by (Wang et al. Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993)

(3.12) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=\unicode[STIX]{x1D705}_{T}\unicode[STIX]{x1D6FB}^{2}T-\frac{{\mathcal{L}}}{c_{p}}\frac{\unicode[STIX]{x2202}p(\unicode[STIX]{x1D719})}{\unicode[STIX]{x2202}t}, & & \displaystyle\end{eqnarray}$$

where the last term corresponds to the release or absorption of latent heat as the phase field varies in time. Note that the fluid is assumed to be at rest in Wang et al. (Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993), but other phase-field models have since included the advection term (Beckermann et al. Reference Beckermann, Diepers, Steinbach, Karma and Tong1999; Anderson, McFadden & Wheeler Reference Anderson, McFadden and Wheeler2000). Using the same non-dimensionalization as in § 2, the phase field (3.6) and the temperature (3.12) read

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D716}^{2}}{m}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}+\frac{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D716}}{St}(\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{M})\frac{\text{d}p}{\text{d}\unicode[STIX]{x1D719}}-\frac{1}{4}\frac{\text{d}g}{\text{d}\unicode[STIX]{x1D719}}, & \displaystyle\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}=-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}-St\frac{\text{d}p}{\text{d}\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}, & \displaystyle\end{eqnarray}$$

where

(3.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}=\frac{{\mathcal{L}}^{2}H\sqrt{a}}{\unicode[STIX]{x1D6FF}c_{p}T_{M}^{2}} & & \displaystyle\end{eqnarray}$$

is the coupling parameter between the phase field and the temperature field. The dimensionless interface thickness and mobility are respectively

(3.16a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}=\frac{\unicode[STIX]{x1D6FF}\sqrt{a}}{H},\quad m=\frac{\unicode[STIX]{x1D6FF}^{2}}{\unicode[STIX]{x1D70F}\unicode[STIX]{x1D705}_{T}}, & & \displaystyle\end{eqnarray}$$

and the Stefan number is defined in (2.8). It is clear from (3.11) that $\unicode[STIX]{x1D716}$ represents the typical interface thickness in the dimensionless space. Since we only consider cases where the bottom boundary is in the liquid phase while the top boundary is in the solid phase, we impose Dirichlet boundary conditions on the phase field

(3.17a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}|_{z=0}=1\quad \text{and}\quad \unicode[STIX]{x1D719}|_{z=1}=0, & & \displaystyle\end{eqnarray}$$

and we recall that we impose the temperature at the boundaries

(3.18a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}|_{z=0}=1\quad \text{and}\quad \unicode[STIX]{x1D703}|_{z=1}=0. & & \displaystyle\end{eqnarray}$$

This phase-field model was initially derived in a much more general context than the classical Stefan problem, focusing on the micro-physics of solidification. It is indeed consistent with the Gibbs–Thompson effects where the temperature at the interface is not exactly the melting temperature, but additionally depends on the local curvature and velocity of the interface. Following the asymptotic analysis of Caginalp (Reference Caginalp1989), Wang et al. (Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993) showed that in the limit of a vanishing interface thickness $\unicode[STIX]{x1D716}\rightarrow 0$ , the following boundary condition applies at the interface

(3.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{M}=-\frac{St}{\unicode[STIX]{x1D6FC}}\left(\unicode[STIX]{x1D705}+\frac{v_{i}}{m}\right), & & \displaystyle\end{eqnarray}$$

where the parameter $St/\unicode[STIX]{x1D6FC}$ can be seen as a dimensionless capillary length, $\unicode[STIX]{x1D705}$ is the dimensionless interfacial curvature and $v_{i}$ is the normal velocity of the interface. Thus, in the additional limit where $St/\unicode[STIX]{x1D6FC}\rightarrow 0$ , and for finite curvature and interface velocity, we recover the original Stefan boundary condition (2.6) where $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{M}$ at the interface, as predicted by Caginalp (Reference Caginalp1989). The value of the mobility is irrelevant provided that the two limits above are respected. In conclusion, the original Stefan problem is recovered provided that

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\ll 1, & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{St}{\unicode[STIX]{x1D6FC}}\ll 1, & \displaystyle\end{eqnarray}$$

while the mobility is fixed to be unity here.

In practice, all the additional parameters introduced by the phase-field formulation are in fact strongly constrained by the limits (3.20)–(3.21). The value of $\unicode[STIX]{x1D716}$ is typically proportional to the numerical grid size in order to accurately solve for the interface region whereas $\unicode[STIX]{x1D6FC}$ is limited by stability constraints. For the interested reader, the effect of these parameters is discussed in more detail in appendix A. In the following, the interface thickness $\unicode[STIX]{x1D716}$ is comparable to the smallest grid size whereas $\unicode[STIX]{x1D6FC}$ is typically of order $St/\unicode[STIX]{x1D716}$ so that both limits (3.20) and (3.21) are satisfied.

3.2 Navier–Stokes equations

The phase-field model described above satisfies the thermal Stefan conditions at the interface given by (2.6)–(2.7). We also have to ensure that the interface corresponds to a no-slip boundary condition for the velocity. Here we choose an immersed boundary method (Mittal & Iaccarino Reference Mittal and Iaccarino2005) called the volume penalization method. The no-slip boundary condition at the liquid–solid interface is implicitly taken into account by adding a volume force to the classical Navier–Stokes equations, solved simultaneously in both liquid and solid domains, leading to

(3.22) $$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D70E}}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}P+Ra\unicode[STIX]{x1D703}\boldsymbol{e}_{z}+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\frac{(1-\unicode[STIX]{x1D719})^{2}\boldsymbol{u}}{\unicode[STIX]{x1D702}}, & & \displaystyle\end{eqnarray}$$

where the last term is the penalization term and $\unicode[STIX]{x1D702}$ is a positive parameter. The incompressibility condition (2.3) is imposed everywhere so that the total volume is necessarily conserved. The penalized (3.22) converges towards the Navier–Stokes equations with a no-slip boundary condition imposed at the interface (Angot, Bruneau & Fabrie Reference Angot, Bruneau and Fabrie1999). The error between the original Navier–Stokes equations and their penalized version scales like $\sqrt{\unicode[STIX]{x1D702}}$ so that $\unicode[STIX]{x1D702}$ is taken as small as possible. This ensures that this term is dominant when $\unicode[STIX]{x1D719}=0$ (i.e. in the solid) and the velocity exponentially decays to zero on a time scale proportional to $\unicode[STIX]{x1D702}$ . When $\unicode[STIX]{x1D719}=1$ (i.e. in the liquid), the penalization term vanishes and the regular Navier–Stokes equations (2.1) are solved. Note that the particular choice $(1-\unicode[STIX]{x1D719})^{2}$ in the numerator of the penalization term is arbitrary and any continuous function that is zero in the liquid and unity in the solid is adequate. For example, in the more general case of porous media, a more complex Carman–Kozeny permeability function $(1-\unicode[STIX]{x1D719})^{2}/\unicode[STIX]{x1D719}^{2}$ can be prescribed, and the momentum and mass conservation equations can also be modified (Beckermann et al. Reference Beckermann, Diepers, Steinbach, Karma and Tong1999; Le Bars & Worster Reference Le Bars and Worster2006). Here we choose the simplest approach of using the phase variable directly to prescribe our penalization term. The quadratic form is chosen in order to match the Carman–Kozeny permeability for $\unicode[STIX]{x1D719}\rightarrow 1$ while recovering a value of unity for $\unicode[STIX]{x1D719}\rightarrow 0$ . Finally, note that all of the results discussed in this paper do not qualitatively depend on this particular prescription of the penalization term; it only affects the detailed structure of the transition between the solid and liquid phase which occurs on length scales typically smaller than the thermal and viscous boundary layers. A comparative study of the different possible expressions for the penalization term, expressed as a function of the phase field itself, is beyond the scope of this paper but would nevertheless prove useful for improving the convergence properties of the current model. The penalization parameter is chosen as small as possible, noting that an explicit treatment of the penalization term leads to stability constraints, typically $\text{d}t<\unicode[STIX]{x1D702}$ (Kolomenskiy & Schneider Reference Kolomenskiy and Schneider2009) where $\text{d}t$ is the time step. In the following, the penalization parameter is chosen so that $\unicode[STIX]{x1D702}=2\text{d}t$ .

We now suppose that the system is two-dimensional which naturally leads to a streamfunction formulation of the Navier–Stokes equations. The streamfunction $\unicode[STIX]{x1D713}$ is defined by $\boldsymbol{u}=-\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D713}\boldsymbol{e}_{y})$ or

(3.23a,b ) $$\begin{eqnarray}\displaystyle u=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\quad \text{and}\quad w=-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}. & & \displaystyle\end{eqnarray}$$

Taking the curl of (2.1) and projecting onto the $y$ -direction leads to the vorticity equation for a two-dimensional flow in the $(x,z)$ plane

(3.24) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}=-\unicode[STIX]{x1D70E}Ra\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}-\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D735}\times (1-\unicode[STIX]{x1D719})^{2}\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{e}_{y}.\qquad & & \displaystyle\end{eqnarray}$$

The no-slip boundary conditions at the top and bottom boundaries correspond to

(3.25a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}|_{z=0,1}=0\quad \text{and}\quad \left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\right|_{z=0,1}=0. & & \displaystyle\end{eqnarray}$$

Similarly, the heat equation, including the phase-field term and the streamfunction decomposition, leads to

(3.26) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}-St\frac{\text{d}p}{\text{d}\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}. & & \displaystyle\end{eqnarray}$$

3.3 Spatial and temporal discretizations

Equations (3.24), (3.26) and (3.13) are solved using a mixed pseudo-spectral finite-difference code. This code has been used in various context from fully compressible convection (Matthews et al. Reference Matthews, Proctor and Weiss1995; Favier & Bushby Reference Favier and Bushby2012) to rapidly rotating Boussinesq convection (Favier, Silvers & Proctor Reference Favier, Silvers and Proctor2014). Each variable is assumed to be periodic in the $x$ -direction and is written as

(3.27) $$\begin{eqnarray}\displaystyle f(x,z)=\mathop{\sum }_{n_{x}}\hat{f}(n_{x},z)\exp (\text{i}k_{x}x)+\text{c.c.}, & & \displaystyle\end{eqnarray}$$

where $n_{x}$ is an integer, $\text{c.c.}$ stands for conjugate terms and the wavenumber is defined as

(3.28) $$\begin{eqnarray}\displaystyle k_{x}=\frac{2\unicode[STIX]{x03C0}n_{x}}{\unicode[STIX]{x1D706}}. & & \displaystyle\end{eqnarray}$$

Horizontal spatial derivatives are computed in spectral space whereas vertical derivatives are discretized using a fourth-order finite-difference scheme.

For the streamfunction, the dissipative fourth-order term is solved implicitly, whereas the advective, temperature and penalization terms are solved explicitly. This is achieved using a classical second-order Crank–Nicolson scheme for the implicit part coupled with a third-order Adams–Bashforth scheme for the explicit part. For the temperature and the phase field, we use a fully explicit third-order Adams–Bashforth scheme. An explicit treatment of these equations is indeed easier due to the nature of the coupling term on the right-hand side of equation (3.26). Note that an implicit scheme could be used to solve these equations (see for example Andersson (Reference Andersson2002)) but the stability constraint associated with solving explicitly both diffusive terms in (3.26)–(3.13) is not very limiting in our two-dimensional case. We have tested the convergence of our numerical scheme in § A.3.

4.1 Critical Rayleigh number

We focus here on the transition between a purely diffusive regime and a convection regime as the fluid depth increases with time. We therefore consider the case where the initial height $h_{0}$ is small enough so that the initial conditions given by (4.1)–(4.3) are stable. It has been shown by Vasil & Proctor (Reference Vasil and Proctor2011) that the convection threshold can be modified compared to the classical Rayleigh–Bénard problem and that a morphological mode grows as soon as $Ra(1-\unicode[STIX]{x1D703}_{M})\overline{h}^{3}>Ra_{c}\approx 1295.78$ . This corresponds to a significant modification of the stability criterion compared to the case of classical no-slip RB convection, for which $Ra_{c}\approx 1707.76$ (Chandrasekhar Reference Chandrasekhar1961). The most unstable wavenumber is also reduced from $k_{c}\approx 3.116$ to $k_{c}\approx 2.552$ . These results are however only valid in the asymptotic limit of large Stefan numbers. While this regime is virtually impossible to reach numerically using the current approach (there exists a time scale separation between the dynamics of the flow and that of the interface), we nevertheless explore this critical transition for finite Stefan numbers in the following.

We start from the initial conditions defined by (4.1)–(4.3) with various initial heights from $h_{0}=0.33$ to $h_{0}=0.45$ , and $\unicode[STIX]{x1D703}_{M}=0.1$ . Using a global Rayleigh number of $Ra=15\,180$ , this leads to an initial effective Rayleigh number varying between $Ra_{e}(t=0)=491$ and $1245$ . We start from infinitesimal temperature perturbations in the liquid layer only. We consider a case where $St=10$ which is large enough to get a reasonable time scale separation while still being accessible numerically. The other numerical parameters are given in table 1 and correspond to case A. We define the kinetic energy density in the system by

(4.9) $$\begin{eqnarray}\displaystyle {\mathcal{K}}(t)=\frac{1}{V_{f}}\int _{V_{f}}\boldsymbol{u}^{2}\,\text{d}V, & & \displaystyle\end{eqnarray}$$

where $V_{f}(t)$ is the volume of fluid as a function of time. The time evolution of the kinetic energy density versus time for various initial heights $h_{0}$ is shown in figure 2. Initially, the kinetic energy in the system briefly increases. This is a consequence of our choice of initial conditions for which the fluid is at rest and temperature perturbations only are added. After this short transient, the kinetic energy density decreases with time for all cases. Surprisingly, the kinetic energy starts to grow for different effective Rayleigh numbers in each case, as early as $Ra_{e}\approx 650$ for the smallest initial height of $h_{0}=0.33$ . The growth rate of this first phase is however much weaker than the typical growth rate when the effective Rayleigh number becomes larger than the classical value of $1707.76$ . We therefore do not observe a clear transition between stable and unstable behaviours at a given critical Rayleigh number, which seems to indicate that perturbations can grow at any value of the effective Rayleigh number.

Figure 2.

Kinetic energy density in the fluid domain versus the time varying Rayleigh number defined by (4.5). Each curve corresponds to a different initial fluid depth $h_{0}$ in (4.2)–(4.3). The dashed curve corresponds to the case where only the horizontally averaged component of the phase field is solved (i.e. the upper boundary is effectively planar).

Table 1.

List of numerical parameters for the different cases discussed in this study.

In order to show that this is a consequence of the upper boundary not being exactly planar, we perform an additional simulation with exactly the same parameters as above and for $h_{0}=0.41$ . The only difference is that we artificially smooth the upper boundary by only solving the horizontally averaged value of the phase field (this is performed in Fourier space by truncating all modes except $k_{x}=0$ ). This is of course artificial but nevertheless useful for understanding the origin of this early growth. The time evolution of the kinetic energy density is shown in figure 2 as a dashed curve. At early times, there are no noticeable differences between the regular simulation and the artificial planar case. However, at later times, there is no growth for the planar case until the critical value of $Ra_{e}\approx 1710$ is reached. This clearly shows that the very early growth of the kinetic energy is associated with the presence of a topography. This topography is very small in amplitude since it is generated by the initial perturbations, but is nevertheless measurable numerically. It is known that any non-planar topography will drive a baroclinic flow at any Rayleigh number (Kelly & Pal Reference Kelly and Pal1978). The amplitude of this gravity current scales linearly with the Rayleigh number and linearly with the amplitude of the topography and is directly forced by the misalignment between the hydrostatic pressure gradient and the inclined temperature gradient normal to the boundary. As we get closer to the threshold $Ra_{e}=1707.76$ , we eventually recover the classical instability mechanism of convection through an imperfect bifurcation (Coullet & Huerre Reference Coullet and Huerre1986). This is consistent with the slow growth of kinetic energy we observed in figure 2, followed by an exponential phase (the growth rate for $Ra_{e}>1707.76$ is actually super-exponential since the Rayleigh number keeps increasing while the instability develops).

There are several reasons why we do not recover the result of Vasil & Proctor (Reference Vasil and Proctor2011) who found a critical transition at $Ra_{e}\approx 1295$ . First, we are not in the asymptotic regime of large Stefan numbers. We repeated the previous simulations at higher $St$ , up to $St=100$ , without qualitative changes in the results discussed above. It is however possible that the asymptotic regime discussed by Vasil & Proctor (Reference Vasil and Proctor2011) is only achieved at much higher Stefan numbers. In addition, Vasil & Proctor (Reference Vasil and Proctor2011) used slightly different boundary conditions and they focused on modes growing on the very slow melting time scale, which is difficult to isolate in our finite Stefan number simulations. Finally, even if we varied all the numerical parameters of our model to confirm that the results discussed above are numerically converged, we cannot discard the possibility that the phase-field approach is inappropriate for studying the evolution of infinitesimal perturbations of the topography, as is the case here. One must remember that the interface is here continuous with a typical width that is here much larger than the perturbations responsible for driving the baroclinic flow. However, we note that once the classical convection instability sets in, all the previous simulations starting from various initial heights lead to the same nonlinear state (apart from a temporal shift as seen in figure 2 at late times) which is discussed in the following sections.

4.2 Nonlinear saturation close to onset and secondary bifurcation

We now explore the evolution of the system once the initial instability saturates and leads to a steady set of convective rolls. In the following, we consider a particular case with a relatively large Stefan number $St=10$ so that we get a reasonable time scale separation between the flow and the interface dynamics. This particular choice is made to simplify the analysis of the bifurcation to be discussed below. For the same reason, we consider a laterally confined case where $\unicode[STIX]{x1D706}=1$ . The other parameters are given in table 1 for case B. We start from an initial height of $h_{0}=0.13$ so that the initial fluid layer is stable ( $Ra_{e}(h_{0})\approx 1245$ ).

After the transient growth discussed in the previous section, we observe at saturation a steady flow and a significant topography which is now clearly non-planar (see figure 3). Perhaps unsurprisingly, the wavelength of this topography is equal to that of the convective rolls below. The solid is locally melting just above rising hot plumes but less so above sinking cold plumes. Once this nonlinearly equilibrated set of rolls and their associated topography exist, the horizontal wavenumber of the rolls is fixed while the average fluid depth keeps increasing with time. This can be seen by measuring the typical horizontal wavelength of the topography as the distance between two local minima of $h(x,t)$ . In order to compare with classical RB convection, we normalize the corresponding wavenumber $k_{x}$ by the time dependent averaged fluid depth $\overline{h}(t)$ . We show in figure 3(a) the effective Rayleigh number of the fluid layer as a function of this normalized wavenumber $\overline{h}k_{x}$ . The marginal stability curve of classical RB with no-slip and fixed temperature walls is shown for reference (Chandrasekhar Reference Chandrasekhar1961). Since the average fluid depth $\overline{h}$ increases, while the horizontal wavenumber of the convection remains constant, the effective Rayleigh number continuously increases like $Ra_{e}\sim \overline{h}^{3}$ . Our simulation closely follows this prediction, as shown in figure 3(a).

Figure 3.

Generation of mean horizontal shear flows during the collapse of steady convective rolls. (a) Effective Rayleigh number as a function of the normalized wavenumber $\overline{h}k_{x}$ . The red curve corresponds to the marginal curve for classical RB (Chandrasekhar Reference Chandrasekhar1961). The oblique dashed lines correspond to a constant horizontal wavenumber and follow $Ra_{e}\sim \overline{h}^{3}$ . (b) Horizontally averaged flow $\overline{u}$ versus $z$ and time. The black line corresponds to the maximum height $\text{max}(h(x,t))$ . The onset of convection and the nonlinear saturation are indicated with vertical dashed lines. (c) Temperature field shown at three successive instants, shown as empty symbols and vertical dotted lines in (a,b). The white line corresponds to the interface defined as $\unicode[STIX]{x1D719}=1/2$ and the grey lines correspond to streamlines.

A simple question now arises: how long can this dynamically evolving set of convective rolls persists against the continuous vertical stretching of the fluid domain? One possibility would be to assume that the rolls are vertically elongated until they become stable again. This is indeed possible since the marginal curve behaves like $Ra_{c}\sim (\overline{h}k_{x})^{4}$ for large wavenumbers whereas the rolls with fixed horizontal wavenumber follow a $Ra_{e}\sim (\overline{h}k_{x})^{3}$ scaling, so that they will eventually become stable as shown in figure 3(a). This is not what is observed in the simulation, however, and a bifurcation occurs well before the possible restabilization of the initially unstable mode. This bifurcation occurs after the rolls have been elongated vertically by an approximate factor of $3$ and corresponds to an abrupt reduction in the horizontal wavenumber $k_{x}(t)$ of the convection rolls.

The detailed nature of this bifurcation can be qualitatively understood by following the time evolution of the horizontally averaged mean flow $\overline{u}(z,t)$ , shown in figure 3(b). This mean flow remains negligible at early times but abruptly grows when the rolls become elongated enough. It first appears as a shear flow with one vertical wavelength, effectively shearing the first set of rolls, as can be seen in the temperature field shown in the middle panel of figure 3(c). Once the initial set of rolls has been disrupted, the mean flow undergoes damped oscillations. A new set of rolls is then generated by the convection instability, with a larger horizontal wavelength, therefore maintaining the unit aspect ratio of the convective cells. This bifurcation is also visible in figure 3(a) where a jump between two $k_{x}\sim \text{const.}$ curves is observed. This is reminiscent of the generation of mean shear flows in laterally confined classical RB convection (Busse Reference Busse1983; Prat, Massaguer & Mercader Reference Prat, Massaguer and Mercader1995; Fitzgerald & Farrell Reference Fitzgerald and Farrell2014; Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014). Our case is however slightly more complicated since the volume of fluid increases with time and the upper topography is non-planar. The generation of mean horizontal shear flows is nevertheless a generic mechanism in RB convection that we also observe in our particular system. Note that the detailed properties of this bifurcation depends on the aspect ratio of numerical domain $\unicode[STIX]{x1D706}$ , the Stefan number $St$ and the nature of the initial perturbations.

This transition between convection rolls of different sizes is expected to repeat itself as the fluid depth keeps increasing, with the additional complication that the effective Rayleigh is ever increasing due to the gradual melting of the solid, so that the bifurcation is expected to become more and more complicated.

4.3 Behaviour at large $Ra$

We now focus on a representative example at high Rayleigh number for which we fix $St=1$ , $Ra=10^{8}$ and $\unicode[STIX]{x1D70E}=1$ . This simulation corresponds to case C in table 1. We consider a domain with a large aspect ratio of $\unicode[STIX]{x1D706}=8$ . Doing so, we aim at minimizing the horizontal confinement effect associated with our periodic boundary conditions. Note that this is the global aspect ratio including the solid domain, the actual aspect ratio of the liquid domain is initially much larger. We therefore expect the liquid phase to display spatio-temporal chaos instead of purely temporal chaos typical of laterally confined systems (Manneville Reference Manneville2006). The initial position of the interface is $h_{0}\approx 0.02$ so that the effective Rayleigh number of the fluid layer is $Ra_{e}(t=0)=Ra(1-\unicode[STIX]{x1D703}_{M})h_{0}^{3}\approx 760$ , well below the critical value. Using a horizontal resolution of $N_{x}=4096$ and a vertical resolution of $N_{z}=1024$ , the smallest grid size is typically $\text{d}x\approx 2\times 10^{-3}$ so that interface width is fixed to $\unicode[STIX]{x1D716}=10^{-3}$ . Assuming a Prandtl–Blasius scaling (Grossmann & Lohse Reference Grossmann and Lohse2000), the dimensionless width of the viscous boundary layers $\unicode[STIX]{x1D6FF}_{v}$ scales like $Ra^{-1/3}$ which typically leads to $\unicode[STIX]{x1D6FF}_{v}\geqslant 2.2\times 10^{-3}$ in our case. Thus, the interface width $\unicode[STIX]{x1D716}$ is always significantly smaller than the viscous boundary layer $\unicode[STIX]{x1D6FF}_{v}$ .

The simulation is run until the interface reaches the upper boundary $z=1$ . For the parameters of this simulation, this approximately takes $0.03$ thermal diffusive time scales. We first show in figure 4 (see also movie 1 in the supplementary materials available online at https://doi.org/10.1017/jfm.2018.773) visualizations showing the temperature and vorticity distributions at different times during the simulation. At the early times, the solution is purely diffusive until the liquid depth reaches its critical value above which convection sets in. Convection is initially steady and laminar, as observed previously, with approximately $132$ convective cells across the whole domain. As the interface progresses, this initial set of convective rolls is vertically stretched, eventually forcing a secondary transition leading to larger convective cells, as discussed previously in § 4.2, although the nature of the bifurcation appears to be different in this large aspect ratio domain (see below). This alternation between quasi-stationary phases of melting where the number of convective cells is conserved and more violent transitions associated with a reordering of the convective cells continues until the upper boundary is eventually reached.

Figure 4.

Visualizations of the total numerical domain for case C in table 1. The temperature is shown on the left (dark red corresponds to $\unicode[STIX]{x1D703}=1$ while dark blue corresponds to $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{M}$ ) while vorticity is shown on the right (blue and red colours correspond to $\pm 0.25\unicode[STIX]{x1D714}_{max}$ respectively). The grey line corresponds to the interface defined by the isosurface $\unicode[STIX]{x1D719}=1/2$ . Time is increasing from top to bottom: $t=5\times 10^{-4}$ , $1.5\times 10^{-3}$ , $6\times 10^{-3}$ , $1.2\times 10^{-2}$ , $2.4\times 10^{-2}$ and $3\times 10^{-2}$ . See also movie 1 in the supplementary materials.

Figure 5.

(a) Position of the interface $h(x,t)$ as a function of $x$ for different times separated by approximately $5\times 10^{-4}$ thermal diffusion times. The colour of the curves corresponds to the signed curvature defined by (4.10). Dark colours correspond to small negative curvatures whereas light colours correspond to cusps with large positive curvatures. (b) Three-dimensional view of the spatio-temporal evolution of the interface $h(x,t)$ . The colour corresponds to the local value of $h(x,t)$ .

Let us first discuss the shape of the interface as the solid continuously melts. We first show in figure 5(a) the interface position as a function of the horizontal coordinate $x$ at different times. The interface is obtained by interpolating the phase field variable in order to find the isocontour $\unicode[STIX]{x1D719}=1/2$ . Equivalently, the interface can be defined as the isotherm $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{M}$ , which leads to the same results, apart from very localized regions of high curvatures where a slight mismatch between the two isocontours is observed, as expected from the Gibbs–Thomson relation (3.19). These discrepancies are however negligible here (i.e. the maximum distance between the isocontours $\unicode[STIX]{x1D719}=1/2$ and $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{M}$ is smaller than the thickness of the boundary layers), as expected from our choice of large coupling parameter $\unicode[STIX]{x1D6FC}$ . The colour of the curves in figure 5(a) corresponds to the signed curvature, derived from the interface position $h(x,t)$ following

(4.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}(x,t)=\frac{\unicode[STIX]{x2202}_{xx}h}{[1+(\unicode[STIX]{x2202}_{x}h)^{2}]^{3/2}}. & & \displaystyle\end{eqnarray}$$

The maximum value of the curvature corresponds to cusps joining two cavities of the topography and is approximately $\unicode[STIX]{x1D705}_{max}\approx 300$ . This is of the same order as the largest curvature achievable by our phase-field approach, which can be approximated by the inverse of the interface width $\unicode[STIX]{x1D716}^{-1}=10^{3}$ . We can therefore be confident that the cusps are numerically resolved and not artificially smoothed by our diffuse interface approach. The horizontal positions of these cusps appear to be very stable, which corresponds to a spatial locking between the convection rolls and the topography (Vasil & Proctor Reference Vasil and Proctor2011). One can also note that these cusps often correspond to small melting rates (i.e. the successive profiles of $h(x,t)$ are close) compared to the much larger cavities with negative curvature where intense localized melting events driven by the underlying hot thermal plumes are observed. An alternative three-dimensional view of the spatio-temporal evolution of the interface is also shown in figure 5(b). The successive bifurcations between different rolls size is clearly visible. We observe various types of cell merging events, from two adjacent cells merging into one, to more complicated behaviours where one cell disappears, leading to the merging of its neighbour cells.

Figure 6.

Spatio-temporal evolution of the temperature at the middle of the fluid layer $\unicode[STIX]{x1D703}(x,\overline{h}/2,t)$ . (a) Full duration of the simulation $0<t<3\times 10^{-2}$ and full spatial extent $0<x<8$ . The dashed line corresponds to the zoom shown below. (b) Zoom in at early times $0<t<6\times 10^{-3}$ ; the dashed lines follow the propagation of a defect in the initially steady array of convective rolls.

We now describe the dynamics of the fluid flow, which is strongly correlated with that of the interface. Our system is not laterally confined so that there is no significant horizontal mean flow, as seen in § 4.2 previously. Instead, we observe a reorganization of the convection cells through local merging events. Figure 6 shows the temperature profile located at the mid-height of the fluid domain, $\unicode[STIX]{x1D703}(x,z=\overline{h}(t)/2,t)$ . At early times, as shown in figure 6(b), the temperature profile is initially purely diffusive and uniform, $\unicode[STIX]{x1D703}(z=\overline{h}/2)=(1-\unicode[STIX]{x1D703}_{M})/2\simeq 1/2$ . The first network of steady convective cells eventually appears and we again observe that the typical horizontal wavelength remains constant after the first nonlinear saturation of the convection instability. The secondary instability, which involved a horizontal mean flow in § 4.2, now appears to be more local since the system is not laterally confined. Interestingly, these local transitions tend to propagate horizontally to neighbouring cells in a percolation process. This is indicated in figure 6(b) by the inclined dashed lines. The typical speed of propagation of this defect in the convection cells lattice, estimated directly from the slope of the dashed lines in figure 6(b), is approximately the same as the fluid vertical velocity. Each cell is therefore destabilized after approximately one turnover time.

In addition, the thermal plumes are clearly oscillatory just after the bifurcation but eventually stabilize and become steady after a short transient. This observation is especially interesting since the Rayleigh number of the system is continuously increasing with time. One might therefore expect the dynamics to become more and more complex, transiting from periodic to chaotic and eventually turbulent solutions as it is the case in classical RB convection without topography. This transition from oscillatory convection to steady convection clearly shows the stabilizing effect that the topography exerts on the flow, locking two counter-rotating convection rolls inside each cavity. At later times and higher Rayleigh numbers, shown in figure 6(a), although the stabilization of the thermal plumes by the topography is still observed, significant temporal fluctuations are nevertheless visible, indicating that the convective cells will eventually transit to more chaotic behaviours. The inevitable transitions between steady, periodic and chaotic solutions observed in classical RB (Gollub & Benson Reference Gollub and Benson1980; Curry et al. Reference Curry, Herring, Loncaric and Orszag1984; Goldhirsch, Pelz & Orszag Reference Goldhirsch, Pelz and Orszag1989) are therefore probably just delayed by the presence of the topography, but will eventually reappear at much larger Rayleigh numbers. This conclusion remains speculative at this stage since this particular simulation is limited to $Ra_{e}<10^{8}$ . It is nevertheless reasonable to expect a different, potentially reduced, interaction between the topography and the underlying flow in the fully developed turbulent regime.

A final interesting observation concerns the clear asymmetry between rising and sinking plumes. Sinking plumes are extremely stable and do not seem to move horizontally, apart from the sudden transitions associated with the reorganization of the convective cells. As seen in figure 4, cold plumes are generated by the merging of two boundary layers descending along the topography, leading to the formation of a high curvature cusp. This cusp is therefore protected by the continuous supply of cold fluid generated by the melting of the neighbouring dome by hot rising fluid. The sinking plumes are therefore always found to be emitted by the cusps. In contrast, rising plumes tend to slowly drift horizontally until they eventually collide with an adjacent sinking plume, leading to destabilization of both convection rolls. The reason for this drift is probably associated with the baroclinic gravity currents, infinitesimal at low Rayleigh numbers and small topography as in § 4.1, but much stronger at large effective Rayleigh numbers and for finite topography slopes. Once a rising thermal plume slightly moves horizontally, it is continuously dragged by the topographic current until a merger occurs. The competition between thermal convection driven by unstable bulk temperature gradients and gravity currents driven by a baroclinic forcing close to an inclined slope is interesting in itself, although we postpone the study of its detailed dynamics to future studies.

5.1 Melting velocity

We now consider the effect of varying the Stefan number on the dynamics of the convective flow and interface. The parameters are the same as previously but we now vary the Stefan number from $St=2\times 10^{-2}$ to $St=50$ . The main consequence of increasing the Stefan number is to increase the time scale separation between the turnover time of the convective cells and the typical time scale of evolution of the topography. As $St$ increases, it takes much more time for a given set of convection rolls to form or alter a topography due to the larger amount of latent heat necessary to do so.

This can be seen in figure 7(a) where the averaged fluid depth is shown versus time for three different Stefan numbers. The dotted lines show the purely diffusive solution in the absence of motions in the fluid phase (i.e. for $Ra_{e}<Ra_{c}$ at all times). These purely diffusive solutions all display the scaling $\overline{h}\sim t^{1/2}$ as expected for diffusive Stefan problems (Vasil & Proctor Reference Vasil and Proctor2011). One observes a departure from this prediction which marks the onset of convection. The larger the Stefan number, the longer it takes to reach the threshold of convection. However, all cases follow the scaling $\overline{h}\sim t$ after the onset of convection. The prefactor however depends on the Stefan number, as shown in figure 7(b). The melting velocity, obtained by a best fit of the previous linear law, is plotted against the Stefan number. The melting velocity seems to be independent of the Stefan number for low values and scales as $St^{-1}$ for large values.

Figure 7.

(a) Time evolution of the horizontally averaged fluid depth $\overline{h}$ for different Stefan numbers. The two scalings $\overline{h}\sim t$ and $\overline{h}\sim t^{1/2}$ are shown as continuous black lines. The two horizontal dotted lines correspond to the critical heights above which convection sets in, estimated from $Ra_{c}\approx 1707$ . (b) Melting velocity $\dot{\overline{h}}$ for different Stefan numbers. The theoretical estimate is derived from (5.5) using $\unicode[STIX]{x1D6FE}\approx 0.115$ and $\unicode[STIX]{x1D6FD}=1/3$ .

These behaviours can be understood by simple energetic arguments. By integrating (3.14) over the whole volume ${\mathcal{V}}$ , we find the relation

(5.1) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\int _{{\mathcal{V}}}[\unicode[STIX]{x1D703}+Stp(\unicode[STIX]{x1D719})]\,\text{d}{\mathcal{V}}=\int _{upper}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}\,\text{d}S-\int _{lower}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}\,\text{d}S, & & \displaystyle\end{eqnarray}$$

where the left-hand side corresponds to the total rate of change of the enthalpy in the system whereas the right-hand side corresponds to the heat fluxes entering and leaving the domain. This conservation constraint must be exactly satisfied at all times during the simulations. Since we work with the temperature as a variable with the latent heat being viewed as an external forcing term, the enthalpy is not explicitly conserved by our scheme. We therefore have to check a posteriori that the enthalpy is indeed conserved in our system. We typically observe a relative error in the total enthalpy of the system of the order of $1\,\%$ at the final time of the simulations when all the solid has melted.

Equation (5.1) can also be used to estimate the rate of change of the average fluid height $\overline{h}(t)$ . The internal heat associated with the solid can be neglected in first approximation since we consider the limit where $\unicode[STIX]{x1D703}_{M}\ll 1$ . In addition, assuming that the fluid layer is fully convective and behaves as in classical RB convection, its average temperature can be approximated by $1/2$ . This leads to the relation (ignoring heat losses from the quasi-isothermal solid above)

(5.2) $$\begin{eqnarray}\displaystyle \left(\frac{1}{2}+St\right)\frac{\text{d}V_{f}}{\text{d}t}=-\int _{lower}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}\,\text{d}S\approx \unicode[STIX]{x1D706}\frac{Nu}{\overline{h}}, & & \displaystyle\end{eqnarray}$$

where the heat flux from the lower boundary has been replaced by the Nusselt number $Nu$ as defined in (4.8) (and we have assumed that $\unicode[STIX]{x1D703}_{M}\ll 1$ ). In (5.2), the volume of fluid

(5.3) $$\begin{eqnarray}\displaystyle V_{f}=\int _{{\mathcal{V}}}p(\unicode[STIX]{x1D719})\,\text{d}{\mathcal{V}} & & \displaystyle\end{eqnarray}$$

can be approximated by $\overline{h}\unicode[STIX]{x1D706}$ and the Nusselt number can be replaced by the usual scaling law involving the effective Rayleigh number of the form $Nu\sim \unicode[STIX]{x1D6FE}Ra_{e}^{\unicode[STIX]{x1D6FD}}$ (see § 5.3 for a more detailed discussion of the heat flux) leading to

(5.4) $$\begin{eqnarray}\displaystyle \left(\frac{1}{2}+St\right)\frac{\text{d}\overline{h}}{\text{d}t}\approx \frac{\unicode[STIX]{x1D6FE}Ra_{e}^{\unicode[STIX]{x1D6FD}}}{\overline{h}}\approx \unicode[STIX]{x1D6FE}Ra^{\unicode[STIX]{x1D6FD}}\overline{h}^{3\unicode[STIX]{x1D6FD}-1}, & & \displaystyle\end{eqnarray}$$

where we have again assumed that $\unicode[STIX]{x1D703}_{M}\ll 1$ . The solution to this equation reads

(5.5) $$\begin{eqnarray}\displaystyle \overline{h}(t)\approx \left[h_{0}^{2-3\unicode[STIX]{x1D6FD}}+\frac{(2-3\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D6FE}Ra^{\unicode[STIX]{x1D6FD}}}{1/2+St}t\right]^{1/(2-3\unicode[STIX]{x1D6FD})}, & & \displaystyle\end{eqnarray}$$

where $h_{0}=h(t=0)$ . Using the typical value of $\unicode[STIX]{x1D6FD}=1/3$ (Grossmann & Lohse Reference Grossmann and Lohse2000), we find that $\overline{h}\sim t$ which is indeed recovered by our simulations. Note that assuming that $\unicode[STIX]{x1D6FD}=1/4$ leads to $\overline{h}\sim t^{4/5}$ which is also in reasonable agreement with our simulations. In addition, assuming that $\unicode[STIX]{x1D6FD}=1/3$ , the melting velocity $\dot{\overline{h}}$ can be estimated for different Stefan numbers from (5.5) and compared with the numerical results. The only adjusting parameter is the prefactor $\unicode[STIX]{x1D6FE}$ linking the Nusselt number with the Rayleigh number. We show in figure 7 the best fit with our numerical data which gives $\unicode[STIX]{x1D6FE}\approx 0.115$ . The agreement is very good over nearly four decades of Stefan numbers. In particular, we recover the fact the melting velocity behaves like $St^{-1}$ in the limit of large Stefan numbers.

5.2 Horizontal and vertical scales of the topography

Let us now discuss some properties of the topography as times evolves. We can track the number of local minima $N_{min}(t)$ of the function $h(x,t)$ as a function of time. The typical length of the cavities generated by flow can be estimated as $l_{c}(t)=\unicode[STIX]{x1D706}/N_{min}(t)$ where $\unicode[STIX]{x1D706}$ is the aspect ratio of the numerical domain. This length scale is plotted in figure 8 for different Stefan numbers as a function of the average fluid depth $\overline{h}(t)$ . As expected, the typical size of the cavities grows with time. Additionally, as seen in figure 4, each cavity contains two convective rolls, each having an opposite circulation. For classical RB convection, the convective rolls typically have a unit aspect ratio (which is related to the fact the critical wavenumber is approximately $k\approx \unicode[STIX]{x03C0}$ ). This corresponds to the typical relation $l_{c}(t)\approx 2\overline{h}(t)$ , which we also plot in figure 8. All the curves remain below this curve, showing that our convective rolls are always slightly elongated in the vertical direction despite the successive dynamical transitions that eventually destabilize them. Note that convective cells tend to be less elongated (i.e. our results get closer to the prediction $l_{c}=2\overline{h}$ ) as the Stefan number increases. This is a direct consequence of the time scale separation typical of the large Stefan number regime, for which the flow can quickly bifurcate to a new, more unstable, set of convective rolls without any significant change in the average fluid depth.

Figure 8.

(a) Typical horizontal extent of the cavities as a function of the averaged fluid depth. The dashed line corresponds to the limit case where each convective roll has a unit aspect ratio. (b) Root-mean-square value of the topography defined by (5.6) as a function of the effective Rayleigh number and for different Stefan numbers. The scaling $Ra_{e}^{1/2}$ is shown for reference.

Finally, it is interesting to consider the typical amplitude of the topography. One can compute the root-mean-square depth as

(5.6) $$\begin{eqnarray}\displaystyle h_{rms}(t)=\sqrt{\overline{(h(x,t)-\overline{h}(t))^{2}}}. & & \displaystyle\end{eqnarray}$$

Figure 8 shows this quantity as a function of the effective Rayleigh number of the fluid layer. Interestingly, there is little dependence with the Stefan number, except close to onset where the saturation of convective occurs later for small Stefan numbers. The typical amplitude of the topography, as measured by $h_{rms}$ , seems to depend mostly on the effective Rayleigh number of the layer, following an approximate scaling of $Ra_{e}^{1/2}$ . Note that the fluctuations of our results, a consequence of the dynamical transitions discussed above, and the relatively small variation in effective Rayleigh numbers, are limiting us from getting a conclusive scaling. It is however interesting to note that the local Reynolds number of RB convection typically scales as $Ra^{1/2}$ (Grossmann & Lohse Reference Grossmann and Lohse2000), so that there might be a link between the amplitude of the topography and the Reynolds number of the underlying flow, and thus the thickness of the boundary layers developing along the topography.

5.3 Effect of the topography on the heat flux

The previous section has shown the non-trivial back-reaction that the topography imprints on the convective flow. The effect of non-uniform boundary conditions on the heat flux in a Rayleigh–Bénard system has a long history. Roughness of the horizontal plates is an obvious candidate to trigger boundary layer instabilities possibly leading to an enhancement of the heat flux (Ciliberto & Laroche Reference Ciliberto and Laroche1999) and possibly to the so-called ultimate regime predicted by Kraichnan (Kraichnan Reference Kraichnan1962; Roche et al. Reference Roche, Gauthier, Kaiser and Salort2010). The wavelength and typical amplitude of our topography is however much larger than the typical roughness used in experiments (Rusaouën et al. Reference Rusaouën, Liot, Castaing, Salort and Chillà2018). Note also that roughness does not always lead to a heat transfer enhancement (Zhang et al. Reference Zhang, Sun, Bao and Zhou2018). In the present case, the horizontal wavelength of the topography is precisely that of the most unstable wavelength of the idealized Rayleigh–Bénard problem in the absence of topography, a situation sometimes referred to as the resonant case (Kelly & Pal Reference Kelly and Pal1978; Bhattacharjee Reference Bhattacharjee1991; Weiss et al. Reference Weiss, Seiden and Bodenschatz2014).

In this section, we consider the heat flux at the bottom boundary $Q_{W}$ defined by (4.7). We show the evolution with time of this heat flux for $St=10$ and $Ra=10^{7}$ as a function of the average fluid depth $\overline{h}$ in figure 9. Before the onset of convection, the purely diffusive heat flux is simply given by $Q_{W}\approx (1-\unicode[STIX]{x1D703}_{M})/\overline{h}$ which is indeed observed initially. After convection sets in, we observe a rapid increase of heat flux associated with the nonlinear overshoot of the instability. The heat flux then tends to decrease with time, but we also observe a succession of plateaus characterized by an approximately constant heat flux $Q_{W}$ , separated by sudden decays. The plateaus correspond to the quasi-steady phases where the convection is locked inside the topography, whereas the sudden decays correspond to the secondary bifurcation where the mean flow is disrupting the convection rolls and inhibits the heat flux across the fluid layer. Starting from the classical relation $Nu\sim \unicode[STIX]{x1D6FE}Ra_{e}^{\unicode[STIX]{x1D6FD}}$ , where $Nu$ is the Nusselt number defined by (4.8), leads to the relation $Q_{W}\sim \overline{h}^{3\unicode[STIX]{x1D6FD}-1}$ so that for $\unicode[STIX]{x1D6FD}<1/3$ , the convective heat flux is indeed a decreasing function of the fluid height, all other parameters being fixed, whereas it is independent of the fluid height when $\unicode[STIX]{x1D6FD}=1/3$ . These different scalings are shown in figure 9 for reference.

Figure 9.

Heat flux averaged over the bottom surface $Q_{W}$ , as defined in (4.7), as a function of the averaged fluid depth $\overline{h}(t)$ . The results are only shown for the case $St=10$ for clarity. The dotted lines correspond to various power law scalings.

We now consider the problematic question of the normalization of this heat flux. In classical RB convection, the diffusive flux across a plane-layer domain is trivially derived from the solution of the purely diffusive heat equation. In our case, however, the diffusive flux is not trivial since the topology of the fluid domain is fully two-dimensional and of finite amplitude. Formally, one should therefore solve the heat equation in order to know the diffusive heat flux across the layer. This refinement has negligible consequence when the topography is of very small amplitude, but this is not the case here, where the topography is of comparable order with the fluid depth. We therefore derive below a second-order correction of the diffusive heat flux at the bottom boundary.

We consider the purely diffusive case of a fluid layer heated from below ( $\unicode[STIX]{x1D703}=1$ in $z=0$ ), with the upper surface at temperature $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{M}$ located at

(5.7) $$\begin{eqnarray}\displaystyle z=h(x,t)=\overline{h}(t)(1+\unicode[STIX]{x1D700}\cos kx), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D700}\ll 1$ , $\overline{h}$ is the mean height given by (4.6), and $k$ is the wavenumber of the topography. We assume that the evolution of $\overline{h}(t)$ is much slower than the diffusion (which is formally justified in the large Stefan number limit: see Vasil & Proctor (Reference Vasil and Proctor2011)). Therefore, we note $\overline{h}=h_{0}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=1-\unicode[STIX]{x1D703}_{M}$ , and we look only for stationary solutions of the diffusion equation

(5.8a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}=0,\quad z\in [0,h(x)], & & \displaystyle\end{eqnarray}$$

with boundary conditions

(5.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}(x,0)=1, & \displaystyle\end{eqnarray}$$

(5.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}(x,h(x))=\unicode[STIX]{x1D703}_{M}. & \displaystyle\end{eqnarray}$$

We expand $\unicode[STIX]{x1D703}$ in power series of $\unicode[STIX]{x1D700}$ ,

(5.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}(x,z)=\unicode[STIX]{x1D703}_{0}(x,z)+\unicode[STIX]{x1D700}\unicode[STIX]{x1D703}_{1}(x,z)+\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D703}_{2}(x,z)+\cdots \,, & & \displaystyle\end{eqnarray}$$

and solve for (5.8) at each order in $\unicode[STIX]{x1D700}$ . After some algebra (cf. appendix B), we obtain at second order

(5.12) $$\begin{eqnarray}\displaystyle Q_{D}=-\frac{1}{\unicode[STIX]{x1D706}}\int _{0}^{\unicode[STIX]{x1D706}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}(x,0)\,\text{d}x=\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D703}}{h_{0}}+\unicode[STIX]{x1D700}^{2}\frac{k\unicode[STIX]{x0394}\unicode[STIX]{x1D703}}{2}\coth kh_{0}. & & \displaystyle\end{eqnarray}$$

On the other hand, the area of the topography per unit horizontal length $A$ (a length per unit length in our two-dimensional geometry) can be computed using (5.7):

(5.13) $$\begin{eqnarray}\displaystyle A=\frac{1}{\unicode[STIX]{x1D706}}\int _{0}^{\unicode[STIX]{x1D706}}\sqrt{1+\left(\frac{\unicode[STIX]{x2202}h(x,t)}{\unicode[STIX]{x2202}x}\right)^{2}}\,\text{d}x=\frac{1}{\unicode[STIX]{x1D706}}\int _{0}^{\unicode[STIX]{x1D706}}\sqrt{1+(-kh_{0}\unicode[STIX]{x1D700}\sin kx)^{2}}\,\text{d}x. & & \displaystyle\end{eqnarray}$$

At order $\unicode[STIX]{x1D700}^{2}$ , we obtain

(5.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}A\equiv A-1={\textstyle \frac{1}{4}}k^{2}h_{0}^{2}\unicode[STIX]{x1D700}^{2}. & & \displaystyle\end{eqnarray}$$

Finally, the diffusion heat flux can be expressed according to this surface area increase

(5.15) $$\begin{eqnarray}\displaystyle Q_{D}=\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D703}}{h_{0}}\left(1+\frac{2\unicode[STIX]{x0394}A}{kh_{0}}\coth kh_{0}\right). & & \displaystyle\end{eqnarray}$$

The typical wavelength in our simulations being of the order of $2h_{0}$ , we can estimate this diffusion flux by writing $k\simeq \unicode[STIX]{x03C0}/h_{0}$ , leading to

(5.16) $$\begin{eqnarray}\displaystyle Q_{D}\simeq \frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D703}}{h_{0}}\left(1+\frac{2\unicode[STIX]{x0394}A}{\unicode[STIX]{x03C0}}\right). & & \displaystyle\end{eqnarray}$$

We can now rescale the heat flux through the fluid layer $Q_{W}$ at each time knowing the averaged fluid depth $\overline{h}$ and the area increase $\unicode[STIX]{x0394}A$ computed numerically at each time step. The diffusive heat flux through the fluid layer is approximately

(5.17) $$\begin{eqnarray}\displaystyle Q_{D}(t)=\frac{1-\unicode[STIX]{x1D703}_{M}}{\overline{h}(t)}\left(1+\frac{2\unicode[STIX]{x0394}A(t)}{\unicode[STIX]{x03C0}}\right) & & \displaystyle\end{eqnarray}$$

and the Nusselt number is finally defined as the ratio between the total heat flux (4.7) and the diffusive flux estimated using (5.17).

Figure 10.

Nusselt number as a function of the effective Rayleigh number $Ra_{e}(t)$ for a Stefan number of $St=10$ . The empty symbols correspond to the classical Rayleigh–Bénard case whereas the dashed line corresponds to the optimum steady solution of Sondak, Smith & Waleffe (Reference Sondak, Smith and Waleffe2015). The thin line is shown for indication and is obtained by computing the effective Rayleigh number on the maximum fluid height instead on its horizontally averaged value.

The result of this normalization is shown in figure 10. For the case $St=10$ , we show the time evolution of the effective Nusselt number as a function of the effective Rayleigh number. For reference, we also show some typical values obtained for classical Rayleigh number (each point corresponding in that case to the time average of a single simulation at fixed Rayleigh number). We also indicate the results of Sondak et al. (Reference Sondak, Smith and Waleffe2015) which correspond to the optimal heat transfer for a two-dimensional steady solution giving $Nu\approx 0.125Ra^{0.31}$ . In Sondak et al. (Reference Sondak, Smith and Waleffe2015), the best fit is actually $Nu\approx 0.115Ra^{0.31}$ but corresponds to a Prandtl number of 7. They obtained a slightly larger prefactor for $\unicode[STIX]{x1D70E}=1$ but the same power law. Interestingly, although our simulation departs significantly from classical Rayleigh–Bénard, our renormalization shows that the Nusselt number follows that of classical RB convection in a quasi-static manner. This is of course only true at large Stefan numbers. For lower Stefan numbers (not shown), the curves are much more erratic and no clear trend can be derived. There are however significant differences between our case and the predictions specific to RB. Although the exponent is not significantly altered by the presence of the melting interface, the prefactor appears larger than what it is for regular RB. This is marginally true at low Rayleigh numbers ( $Ra_{e}<10^{6}$ ) but quite clear at higher Rayleigh numbers. This can be attributed to the back reaction of the topography on the convective rolls, which appears to be a stabilizing effect by delaying the transition to periodic convection. In two dimensions, this transition typically occurs around $Ra\approx 10^{5}$ and reduces the Nusselt number, as can be seen in figure 10. The presence of the topography induced by the convective flow itself seems to favour stable quasi-steady rolls as opposed to oscillatory ones. This leads to an increase in heat flux when compared to classical RB and is closer to the optimal solution of Sondak et al. (Reference Sondak, Smith and Waleffe2015), derived assuming steady laminar solutions. This marginal increase in the Nusselt number was also recently reported in the independent study by Rabbanipour Esfahani et al. (Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018), both in two and three dimensions, although for a Prandtl number of $10$ . Note finally that although we carefully normalized the heat flux, our choice of Rayleigh number is rather arbitrary. One could argue for example that the effective Rayleigh number based on the averaged fluid depth is barely relevant and that only the maximum depth where the Rayleigh number is effectively maximum matter for the heat flux. The results corresponding to this particular choice are shown in figure 10 as the thin line. Although it does reduce the overall heat flux for a given Rayleigh number, our conclusions drawn above remain qualitatively valid. Since we are limited in the maximum value for our effective Rayleigh number (a consequence of the finite vertical extent of our numerical domain), it is not clear how the topography affects the heat flux in the fully chaotic regime reached at much higher Rayleigh numbers.