Melting driven by rotating Rayleigh–Bénard convection

Abstract

We study numerically the melting of a horizontal layer of a pure solid above a convecting layer of its fluid rotating about the vertical axis. In the rotating regime studied here, with Rayleigh numbers of order $10^{7}$, convection takes the form of columnar vortices, the number and size of which depend upon the Ekman and Prandtl numbers, as well as the geometry – periodic or confined. As the Ekman and Rayleigh numbers vary, the number and average area of vortices vary in inverse proportion, becoming thinner and more numerous with decreasing Ekman number. The vortices transport heat to the phase boundary, thereby controlling its morphology characterized by the number and size of the voids formed in the solid, and the overall melt rate, which increases when the lower boundary is governed by a no-slip rather than a stress-free velocity boundary condition. Moreover, the number and size of voids formed are relatively insensitive to the Stefan number, here inversely proportional to the latent heat of fusion. For small values of the Stefan number, the convection in the fluid reaches a slowly evolving geostrophic state wherein columnar vortices transport nearly all the heat from the lower boundary to melt the solid at an approximately constant rate. In this quasi-steady state, we find that the Nusselt number, characterizing the heat flux, co-varies with the interfacial roughness, for all the flow parameters and Stefan numbers considered here. This confluence of processes should influence the treatment of moving boundary problems, particularly those in astrophysical and geophysical problems where rotational effects are important.

1. Introduction

The coupling between a solid and the liquid from which it forms controls the long term fate of both phases. Through deliberate manipulation of the flow of the nutrient phase, engineers aim to control the character of a solidified material (Davis 2001). When the heat transport required for solidification occurs through diffusion, initially planar phase boundaries remain planar. But the presence of convection invariably leads to non-planar interfaces. The uncontrolled interplay of convection, rotation and phase change determines the dynamics of many geophysical and astrophysical systems. Indeed, such processes operate from the Earth's core to the principal components of the cryosphere (e.g. Huppert 1990; Worster 2000). In astrophysics, they underlie planet formation (e.g. Armitage 2020), wherein for example the proto-Earth was believed to rotate approximately ten times faster than today (e.g. Cuk & Stewart 2012), and the growth of neutron star crusts (e.g. Baym et al. 2018), amongst many other phenomena. The confluence of dynamic and thermodynamic processes in such systems is highly complex and involves multiple time scales, components and phases.

Here, we study a simplified system of a single-component rotating phase boundary heated from below. The associated rotation-influenced convection brings heat to the solid upper boundary, controlling the morphology of the melting solid.

A non-rotating layer of fluid heated from below begins convecting when the thermal buoyancy overcomes the viscous and thermal dissipation effects that suppress vertical motions. This balance is characterized by the dimensionless Rayleigh number

(1.1)\begin{equation} Ra=\frac{g\alpha\Delta Th^{3}}{\nu\kappa}, \end{equation}

where $g$ is the acceleration due to gravity; $\alpha$, $\nu$ and $\kappa$ are the coefficient of thermal expansion, the viscosity and the thermal diffusivity of the fluid; and $h$ is the depth of the fluid layer across which a temperature difference $\Delta T$ is imposed. Convective motions begin when $Ra$ exceeds a critical value $Ra_{c} = {O}(10^{3})$, the prefactor depending on the boundary conditions (e.g. Chandrasekhar 1961).

In direct analogy with stratification in non-rotating systems, rotation suppresses vertical motions due to buoyancy (Veronis 1970). Therefore, the critical Rayleigh number above which convection occurs is a function of the rotation rate of the system (Chandrasekhar 1953). The Ekman number is the relevant non-dimensional rotation rate and is

(1.2)\begin{equation} E=\frac{\nu}{2\varOmega h^{2}}, \end{equation}

where $\varOmega$ is the angular velocity of the system. Thus, rapidly rotating systems are characterized by small $E$. Whereas in non-rotating convection a given set of boundary conditions determines the single value of $Ra_{c}$, in rotating convection $Ra$ is an increasing function of $E^{-1}$, where both the functional form and numerical factors depend on the boundary conditions of the problem.

If the horizontal directions are assumed to be periodic, the onset of convection occurs above $Ra_{c}\sim E^{-4/3}$. For one free-slip and one no-slip boundary each (and periodic boundary conditions in the horizontal), in the limit of large $E^{-1}$ (Chandrasekhar 1953), $Ra_{c}$ is

(1.3)\begin{equation} Ra_{c}^{{bulk}}=2.39E^{{-}4/3}. \end{equation}

If the horizontal directions are bounded by walls, the critical Rayleigh number for the so-called ‘wall mode’ (Zhong, Ecke & Steinberg 1991; Ecke, Zhong & Knobloch 1992) is, in the limit of large $E^{-1}$, given by (Herrmann & Busse 1993)

(1.4)\begin{equation} Ra_{c}^{{wall}}|_{E^{{-}1}\rightarrow\infty}={\rm \pi}^{2}(6\sqrt{3})^{1/2}E^{{-}1} < Ra_{c}^{{bulk}}. \end{equation}

In a rotating system bounded laterally by walls, flow is absent for $Ra < Ra_{c}^{{wall}}$. The flow structures that appear for $Ra>Ra_{c}^{{wall}}$ take the form of a peripheral streaming current adjacent to the walls, with alternating bands of up- and down-welling flow. While the flow in them is still cyclonic, these patterns precess about the axis of rotation in a retrograde direction (Horn & Schmid 2017; De Wit et al. 2020; Favier & Knobloch 2020; Zhang et al. 2020), even when there are severe obstacles in the way (Favier & Knobloch 2020). These wall modes persist even when the bulk of the flow begins to convect, and they underlie an observed mismatch between theoretical and numerical predictions of heat transport and laboratory observations at large $Ra$ (De Wit et al. 2020).

When $Ra>Ra_c^{{bulk}}$, convection begins throughout the fluid. For $Ra\lesssim 10 Ra_{c}^{{bulk}}$, flow occurs along columnar (Taylor) vortices that span the depth of the fluid (Boubnov & Golitsyn 1986, 1990; Zhong et al. 1991; King et al. 2009; Aurnou et al. 2015). These vortices are predominantly cyclonic near the upper and lower boundaries, with equal numbers of cyclonic and anticyclonic vortices in the interior (Boubnov & Golitsyn 1986; Zhong et al. 1991; Vorobieff & Ecke 1998; Kunnen, Clercx & Geurts 2010), thereby transporting heat from the boundaries (Sakai 1997). For $Ra>10 Ra_{c}^{{bulk}},$ the columnar vortices become plume-like and lose their vertical alignment with the axis of rotation. The highest Rayleigh numbers achieved in our simulations are in this regime. For sufficiently large $Ra$ (and sufficiently large $E^{-1}$), a state of ‘geostrophic turbulence’ sets in (Boubnov & Golitsyn 1990; King, Stellmach & Aurnou 2012; Shi et al. 2020), a computationally challenging regime to study.

The nature of rotating convection and the rate of heat transport are controlled by the combination of $E$, $Ra$ and $Pr$, and thus so too will be the melt rate and patterns of an adjacent phase boundary, such as we study here. While varying the dimensionless latent heat, or Stefan number, is expected to affect the overall rate of phase change, the effects on the interfacial patterns that form are more subtle, which largely reflect the nature of the transport properties of the bulk flow. This confluence of effects forms the core of our study.

The rest of the paper is organized as follows. We describe the structure of the problem in § 2, providing details of the phase-change treatment used; the approximations made; the relevant physical scales and the non-dimensionalization; the boundary and initial conditions; and the numerical algorithm used to solve the governing equations. In § 3, we discuss the effects of the control parameters on the phase boundary morphology, which is dominated by rotation. We obtain the melt rates and their associated Nusselt number dependencies. Additionally, we discuss how the dynamics changes if the system is periodic in the horizontal, if the lower boundary is one of no slip and when the solid has a thermal diffusivity different from the liquid. We conclude with some ideas for future work.

2. Structure of the problem

Our study geometry is a box of dimensions $L\times L\times H$, with gravity $g$ in the $-z$ direction, and rotating about the $+z$ axis with an angular velocity $\varOmega$, shown schematically in figure 1. The aspect ratio of the simulation domain is $L/H=2$. The mean height of the liquid layer at time $t$ is $h(t)$, with $h(t=0)=h_{0}$. We use the domain half-height as our length scale (see § 2.2 below), and define the aspect ratio as $A=2L/H$. The system is heated from below by imposing a constant temperature difference between the lower and upper boundaries, thereby melting the solid. As described in § 2.3, the majority of our results are obtained with the entire solid at the melting temperature, so that there is no heat conduction through the solid.

Figure 1. (a) A schematic of the geometry used, with the coordinate directions and dimensions marked. The initial liquid height is $h(t=0) = h_0$. (b) Vertical cross-section of the geometry considered at $t>0$. The system rotates about the $z$ axis, and gravity is in the $-z$ direction. Here, $T_m$ is the melting temperature of the pure substance, and the lower boundary is at temperature $T_m + \Delta T$. The effective Rayleigh and Ekman numbers are defined based on the horizontally averaged fluid height $h(t)$, while the reference values are defined based on $H/2$, where $H$ is the height of the solid+liquid system.

2.1. Enthalpy method

We employ a mixture theory approach to tracking the solid region, such that a solid fraction variable $\chi$ varies from $0$ in the liquid state to $1$ in the solid state. The densities of the solid and liquid phases are $\rho _{s}$ and $\rho _{l}$ respectively; their heat capacities are $C_{s}$ and $C_{l}$ respectively; and the latent heat of fusion is $\lambda$. Here, for simplicity, we only consider the case where the solid and liquid have the same densities and

(2.1)\begin{gather} \rho_{s} =\rho_{l}\left(=\rho\right), \end{gather}

(2.2)\begin{gather}C_{s} =C_{l}\left({=}C_{p}\right), \end{gather}

with $\rho$ and $C_{p}$ being constants. The solid and liquid enthalpies are

(2.3)\begin{gather} \mathcal{H}_{s} =\rho C_{p}T,\text{ and } \end{gather}

(2.4)\begin{gather}\mathcal{H}_{l} =\rho C_{p}T+\rho\lambda, \end{gather}

respectively. The enthalpy of the solid phase at the melting temperature $T_{m}$ is $\mathcal {H}_{0}=\rho C_{p}T_{m}$, and that of a mixture of solid and liquid phases with solid volume fraction $\chi$ is given by

(2.5)\begin{align} \mathcal{H} & =\chi\rho C_{p}T+\left(1-\chi\right)\rho\left[C_{p}T+\lambda\right].\nonumber\\ & =\rho C_{p}T+\left(1-\chi\right)\rho\lambda . \end{align}

We non-dimensionalize the enthalpy as

(2.6)\begin{equation} \phi=\frac{\mathcal{H}-\mathcal{H}_{0}}{\rho C_{p}\Delta T}=\frac{T-T_{m}}{\Delta T}+\frac{\lambda}{C_{p}\Delta T}\left(1-\chi\right), \end{equation}

where $\Delta T$ is the difference between the temperature of the lower boundary and the melting temperature. Thus, if

(2.7)\begin{equation} \theta = \frac{T-T_{m}}{\Delta T} \end{equation}

is defined to be the non-dimensional temperature, and

(2.8)\begin{equation} St=C_{p}\Delta T/\lambda\end{equation}

is the Stefan number (often also defined as the inverse of this) then we have

(2.9)\begin{equation} \phi=\theta+St^{{-}1}\left(1-\chi\right). \end{equation}

We note that in the purely solid state $\chi =1$ and ${\theta \leq 0},$ so that $\phi \leq 0$. The equation of state (2.9) can be inverted to give the solid fraction in terms of the enthalpy as

(2.10)\begin{equation} \chi=1-\max\left[0,\min\left(1,St\, \phi\right)\right], \end{equation}

and hence the temperature follows as

(2.11)\begin{equation} \theta=\phi-St^{{-}1}\left(1-\chi\right). \end{equation}

Thus, in a pure solid, $\chi =1$, $\theta =\phi$; in a pure liquid, $\chi =0$, $\theta =\phi -St^{-1}$; in the mixed phase, $0<\chi <1$ and $\theta = 0$, by definition. In the vicinity of the phase boundary $\chi$ must change from 0 to 1 over a very thin region (see e.g. Rabbanipour Esfahani et al. 2018; Favier, Purseed & Duchemin 2019), which is a requirement that our simulations obey. The normal motion of the phase boundary, $u_{m}$, is determined by the interphase difference between heat fluxes, and the Stefan condition in dimensional variables is

(2.12)\begin{equation} \rho\lambda u_{m}=k_{s}\left(\boldsymbol{\nabla} T\right)_{s}-k_{l}\left(\boldsymbol{\nabla} T\right)_{l}, \end{equation}

where $(\boldsymbol {\nabla } T)_{s}$ and $(\boldsymbol {\nabla } T)_{l}$ are the temperature gradients in the solid and the liquid on either side of the phase boundary; and $k_{s}$ and $k_{l}$ are the thermal conductivities in the solid and liquid respectively.

2.2. Governing equations

The equations of motion that govern the evolution of the velocity $\boldsymbol {u}$, and the enthalpy $\phi$, defined in (2.9), are as follows. We study the rotating Oberbeck–Boussinesq equations with the assumptions in (2.1) and (2.2), which are

(2.13)\begin{gather} \frac{\textrm{D}\boldsymbol{u}}{\textrm{D}t} ={-}\frac{\boldsymbol{\nabla} p}{\rho}+\nu\nabla^{2}\boldsymbol{u}+g\alpha\boldsymbol{e}_{z}\left(T-T_{m}\right)-2\varOmega\boldsymbol{e}_{z}\times\boldsymbol{u}, \end{gather}

(2.14)\begin{gather}\frac{\textrm{D}\theta}{\textrm{D}t} = \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\kappa\boldsymbol{\nabla}\theta\right), \text{ and } \end{gather}

(2.15)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}

where ${\rm D}/{\rm D}t$ is the material derivative, $\alpha$ is the coefficient of thermal expansion, $\nu$ is the kinematic viscosity of the fluid and $\kappa =\chi \kappa _{s}+(1-\chi )\kappa _{l}$ is the local thermal diffusivity. These equations are non-dimensionalized using the temperature scale $\Delta T$ from (2.9), and the length scale $H/2$, where $H$ is the height of the domain, giving a buoyancy velocity $U_{b}=(g\alpha \Delta T H/2)^{1/2}$. Using these scales, the dimensionless equations of motion become

(2.16)\begin{gather} \frac{\textrm{D}\boldsymbol{u}}{\textrm{D}t} ={-}\boldsymbol{\nabla} p+\left(\frac{Pr}{Ra}\right)^{1/2}\nabla^{2}\boldsymbol{u}+\boldsymbol{e}_{z}\theta-{Ro_c^{{-}1}}\boldsymbol{e}_{z}\times\boldsymbol{u}, \end{gather}

(2.17)\begin{gather}\frac{\textrm{D}\theta}{\textrm{D}t} = \left(\frac{1}{Ra Pr}\right)^{1/2}\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\hat{\kappa}\boldsymbol{\nabla}\theta\right),\text{ and} \end{gather}

(2.18)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}

where $Pr=\nu /\kappa _{l}$ is the Prandtl number, $Ro_c$ is the Rossby number (see (2.22)) and $\hat {\kappa }=\kappa /\kappa _{l}$ is the ratio of the local thermal diffusivity to the diffusivity in the liquid. The Stefan condition ((2.12)) in non-dimensionalized form is given by

(2.19)\begin{equation} u_{m}=\frac{St}{Re\cdot Pr}\left[\hat{\kappa}_{s}\left(\boldsymbol{\nabla}\theta\right)_{s} -\left(\boldsymbol{\nabla}\theta\right)_{l}\right],\end{equation}

where $\hat {\kappa }_{s}=\kappa _{s}/\kappa _{l}$ is the non-dimensional thermal diffusivity in the solid. Finally, in the solid there is only heat conduction and hence $\boldsymbol {u}=0$ in ((2.16) and (2.17)).

As the solid melts and the height of the liquid layer increases, the effective Rayleigh and Ekman numbers evolve according to

(2.20)\begin{gather} Ra_{{eff}} = Ra\left[\frac{h\left(t\right)}{H/2}\right]^{3},\quad\text{and} \end{gather}

(2.21)\begin{gather}E_{{eff}} = E\left[\frac{H/2}{h\left(t\right)}\right]^{2}, \end{gather}

respectively, showing that as the solid melts and the liquid layer becomes deeper, $Ra_{ {eff}}$ and $E^{-1}_{{eff}}$ both increase. We also note that the ratio $(Ra/Ra_{c}^{{bulk}})_{{eff}}\sim RaE^{4/3}$ (from (1.3)) increases with time as $h^{1/3}$.

Unless specifically mentioned, we label the results presented here with the reference values $Ra$ and $E$. The effective Rayleigh and Ekman numbers $Ra_{{eff}}$ and $E_{{eff}}$ are considered in the heat transport calculations in § 3.4.

Lastly, the Rossby number $Ro_c$ in (2.16), also sometimes called the convective Rossby number, is a measure of the rotation dominance of the flow, and is given by

(2.22)\begin{equation} Ro_c = \left(\frac{Ra}{Pr Ta}\right)^{1/2}= E\left(\frac{Ra}{Pr}\right)^{1/2}, \end{equation}

where $Ta=E^{-2}$ is the Taylor number. Despite system specific definitions of the Rossby number, such as in geophysical fluid dynamics (see e.g. Cushman-Roisin & Beckers 2011, chapter 9), all flows with Rossby numbers much less than unity are rotationally dominated.

2.3. Initial and boundary conditions

At $t=0$, both the solid and liquid phases are at the melting temperature $\theta =0$. Unless otherwise mentioned, we use $h_{0}=H/2$. The upper and lower boundaries are held at temperatures $\theta =-f$ and $\theta =1$ respectively ($f=0$ except in § 3.2.3). The lateral boundaries are insulating, no-slip walls. No-slip conditions are also applied at the freely evolving phase boundary, where the temperature is $\theta =0$. Ravichandran & Wettlaufer (2020) showed that free-slip boundaries support flow structures that no-slip boundaries cannot. Here, in order to examine how such structures influence the melting dynamics, a free-slip velocity condition is used on the lower boundary, except in § 3.2.1, where we study the influence of the no-slip velocity boundary condition on the lower boundary.

2.4. Numerical simulations

Equations ((2.16) and (2.17)), together with (2.19), are solved using the finite volume solver Megha-5 on a uniform grid in all three space directions (Diwan et al. 2014; Prasanth 2014; Ravichandran & Wettlaufer 2020; Ravichandran, Meiburg & Govindarajan 2020). After every time step of (2.16) and (2.17), an equilibration step is implemented using (2.10) and (2.11). This procedure is similar to that used by Rabbanipour Esfahani et al. (2018) and has been validated against analytical results (appendix A). The requisite velocity conditions in the resulting arbitrarily shaped solid region are applied using the volume penalization method of Kevlahan & Ghidaglia (2001), wherein the solid is modelled as a porous medium with vanishing porosity. This amounts to adding a term $-({\chi }/{\eta })\boldsymbol {u}$ to the right-hand side of (2.16), where $\eta \ll 1$ is the penalization parameter. Our simulations are performed with up to $512^{2}\times 256$ grid points, a penalization parameter of $\eta =2\times 10^{-3}$, and a time step of $\delta t = 10^{-3}$. The results presented are independent of the grid resolution and insensitive to the value of the penalization parameter used (appendix B).

We note that for the single component two-phase system considered here, the solid–liquid interface has to be sharp and hence $\chi$ varies smoothly from $0$ to $1$ over a finite number of grid points (see figure 23 in appendix A). For the purposes of plotting, the solid–liquid interface is taken to be the iso-surface $\chi =0.5$.

3. Results and discussion

The range of Ekman and Rayleigh numbers we consider here are listed in table 1, and correspond to rapidly rotating convection. For the associated values of $Ra/Ra_{c}$, we obtain no flow for $Ra < Ra_{c}^{{wall}}$; a streaming flow close to the walls (the ‘wall modes’) for $Ra_{c}^{{wall}} < Ra < Ra_{c}^{{bulk}}$; and columnar vortices for $Ra>Ra_{c}^{{bulk}}$. We do not study the geostrophic turbulence regime, $Ra\gg Ra_{c}^{{bulk}}$. In the majority of cases we report here, the flow takes the form of columnar vortices, with a peripheral retrograde near-wall current. We show how the nature of the flow controls the morphology of the melting of the solid, and how the melting influences the flow structures. We also study the sensitivity of these results to the Stefan number. As we explain below, choosing a Prandtl number of $5$ allows columnar vortices to form at lower Rayleigh numbers.

Table 1. (a) Ranges of the controlling parameters (defined in the text) used. (b) Typical boundary conditions or values of parameters used, except in special cases called out in the text.

3.1. Flow structure and melting morphology

Before discussing the influence of these flow structures on the morphology of the melting, we look first at some properties of the columnar vortices themselves. We identify these vortices as isolated regions at the horizontal plane given by $z=H/4$ where

(3.1)\begin{equation} \omega_{z}=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}> \omega_{0}. \end{equation}

Whilst the threshold used, $\omega _{0}=0.25$, is arbitrary, this choice does not change the number of vortices significantly, but it does affect the vortex area, as is to be expected. The rotating convection driving the melting is time dependent, and the mean and maximum vorticity increase with time. For this reason, we rationalize an arbitrary threshold in order to have a means of comparing vortex areas and numbers at different points of time.

3.1.1. Rotational dominance and columnar vortices

For a given $Ra, Pr$ combination, decreasing $E$ increases the rotational control of the flow and we expect a larger number of thinner vortices (Zhong et al. 1991; Sakai 1997; Vorobieff & Ecke 1998), as shown in figure 2(a). Moreover, as the Rayleigh number increases the number of vortices decreases, as shown in figure 2(b). Of particular relevance to the phase-change dynamics, figure 3(a) shows that as the number of vortices increases the average area of each vortex decreases. Moreover, this behaviour is independent of the flow regimes studied, as evidenced by the parametric collapse onto a single curve. Figure 3(b) shows that, beyond the initial transients, the total vortex area reaches a quasi-steady state. For a given $E$, this total vortex area increases with increasing $Ra$.

Figure 2. The number of columnar vortices as a function of time for $Pr=5$, $St=1$ and (a) $Ra=7.8\times 10^{6}$, and (b) $E=10^{-4}$, showing that, as rotational effects become more dominant, the number of vortices increases.

Figure 3. Dependence of flow structure on the flow parameters $E$ and $Ra$ with $Pr=5$ and $St=1$. (a) The number of vortices and the average area of each vortex area inversely proportional to each other. (b) The total cross-sectional area of the columnar vortices is an increasing function of time before saturating at late times.

Vertical and horizontal cross-sections of the temperature and vertical velocity in figure 4 show the typical patterns of flow and melting seen at the smallest and largest $E$ in our simulations (table 1). Particularly notable is the increase in the number of vortices in smaller $E$ more rotationally dominant flows.

Figure 4. Cross-sections of the temperature $\theta$ and the vertical velocity $w$ for (a–d) $E=10^{-3}$, $Ra=2\times 10^{5}$, $Pr=5$, $f=0$, $St=1$, $t=240$; and (e–h) $E=8\times 10^{-5}$, $Ra=7.8\times 10^{6}$, $Pr=5$, $f=0$, $St=1$, $t=500$. In each panel, the horizontal sections (a,b,e,f) are plotted on the $z=H/4$ plane and the vertical sections (c,d,g,h) are plotted on the $y=0$ plane. The yellow lines in the vertical sections show the instantaneous location of the solid–liquid interface. Vertical heat transport occurs in columnar vortices as reflected in the pattern of the melting solid.

These columnar vortices carry heat from the lower boundary to the solid and, as figure 4 shows, etch voids into the solid. Therefore, the morphology of the phase boundary – the average area and number of void regions melted into the solid – reflects the state of the flow. Figure 5 shows that the number of voids and their average cross-sectional area are proportional to the number and the average area of the vortices respectively. However, whereas the number and size of the vortices play a role in the total heat transport by the fluid, the heat transfer is not simply proportional to the total vortex area, but depends additionally upon their specific heat and velocity, as described presently.

Figure 5. (a) The number of solid voids as a function of the number of vortices, showing the linear dependence of the former on the latter. (b) The area of the solid voids as a function of the area of the vortices. In both figures, points are plotted every $20$ flow units excluding initial transients and before the fluid comes into direct contact with the upper boundary. Here, $Pr=5$, $f=0$ and $St=1$ in all cases shown.

We note that figure 4 shows sharp cusps in the solid–liquid interface. Such cusps are a common challenge in numerical simulations of interfacial flows (e.g. Popinet 2018). Here, we find no evidence that these features influence the overall dynamics appreciably. In particular, we have verified that the shapes and sizes of the cusps, and the shapes and areas of the voids are independent of grid resolution.

As the melting proceeds and the height of the liquid layer grows, $Ra_{{eff}}$ and $E_{{eff}}^{-1}$ grow as well ((2.20) and (2.21)). Moreover, as the vortices merge into larger vortices, the voids do as well. The average area of the voids thus grows as a function of time, as seen in the plot of the average void area vs $E_{{eff}}$ in figure 6. We note, however, that figure 6 is primarily intended to motivate future work. Namely, because they do not span two decades on both axes, a rigorous evaluation (see, e.g. Stumpf & Porter 2012) of the relationship between the void area and $E_{{eff}}$ cannot be made.

Figure 6. (a) The average area of the voids formed grows with time, and is seen to grow proportionally to $E_{{eff}}$. Apart from the initial transients (and the divergence to infinity in cases where all the solid has melted away within the simulation time), the same proportionality holds for different values of $E={E_{{eff}}(t=0)}$. (b) A reasonable collapse is obtained if the void areas are multiplied by $E^{-3/2}$ (note that we multiply by the initial value, not the abscissa). The parameter combinations are the same as in figure 5, and the symbols have the same meaning.

The convective Rossby number, $Ro_{c}$, is another key parameter that quantifies the rotational control of the flow. As seen in (2.22), for a given combination of $Ra$ and $E$, a larger $Pr$ leads to a smaller $Ro_{c}$, and thus to greater rotational dominance. In figure 7, this is reflected in the melt voids that are created by the columnar vortices present for $Pr=5$, but absent for $Pr=1$.

Figure 7. The solid–liquid interface (viewed from the solid side) with $St=1$, $f=0$ for (a) $E=10^{-3}$, $Ra=2\times 10^{5}$, $Pr=1$, $t=120$; (b) $E=10^{-3}$, $Ra=2\times 10^{5}$, $Pr=5$, $t=240$; (c) $E=10^{-4}$, $Ra=5\times 10^{6}$, $Pr=1$; $t=240$; (d) $E=10^{-4}$, $Ra=5\times 10^{6}$, $Pr=5$, $t=500$. For $Pr=5$, vertical heat transport occurs in columnar vortices as reflected in the pattern of the melting solid.

We note that the times at which the phase boundaries are shown in figure 7 reflect that for a given $Ra$, a reduction in $Pr$ reflects an increase in heat transfer and hence melt rate, further evidence of which is seen in figure 8, where we plot the amount of solid $h_{s}(t)=H-h$ as a function of time for $Pr=1$ and $Pr=5$.

Figure 8. The volume-averaged height of the solid $h_{s}=H-h$ as a function of time, showing the role of the flow parameters, with $St=1$, $f=0$. (a) $Ra = 10^5;$ Columnar vortices are absent for both Prandtl numbers. (b) $Ra = 10^6;$ Columnar vortices are present for $Pr=5$. For the three combinations of $E,Ra$ (i) $E=10^{-3}$, $Ra=10^{5}$, $Ra/Ra_{c}^{{bulk}}=4.2$; (ii) $E=5\times 10^{-4}$, $Ra=10^{5}$, $Ra/Ra_{c}^{{bulk}}=1.6$; (iii) $E=5\times 10^{-4}$, $Ra=10^{6}$, $Ra/Ra_{c}^{{bulk}}=16.6$. For a given $Ra$, melting is slower for larger $Pr$ regardless of the degree of supercriticality $Ra/Ra_{c}$ or the presence of columnar vortices.

It is intuitive that for a given $E$, the melt rate increases with $Ra$ and this is seen in figure 9(a,b). Moreover, for similar values of $Ra/Ra_{c}$, melting is faster for larger $E$, when vertical transport is less rotationally constrained. We analyse the energy balance underlying the melting rates and the effective Nusselt numbers in detail in § 3.4.

Figure 9. The volume-averaged height of the solid $h_{s}$ as a function of time, showing the role of the Ekman and Rayleigh numbers for $Pr=5$, $f=0$ and $St=1$ (a–c) and $St=0.2$ (d). Increasing $Ra/Ra_{c}$ leads to a larger melt rate, as shown for (a) $E=10^{-3}$ and (b) $E=5\times 10^{-4}$. For comparable $Ra/Ra_{c}$, melting is slower for smaller $E$, as seen in (c,d). Note that the simulations in (d) are run for $2000$ flow time units.

3.1.2. Wall modes and peripheral melting

When the Rayleigh number approaches the critical value, $Ra_{c}^{{bulk}}$, heat is transported predominantly through the peripheral streaming current, and hence the solid regions closer to the walls melt significantly faster than the interior, which, as shown in figure 10, remains more planar. Whereas in figure 10(a), $Ra/Ra_{c}$ = ${O}(1)$, as it increases we see both the effects of the wall modes and the bulk flow. Thus, when columnar vortices are present, as is the case for $Pr=5$ in figure 10(b), the voids formed penetrate deeper into the solid than the melt regions created by the wall modes.

Figure 10. The height of the liquid layer, averaged in the horizontal direction for $y\in [-H/2,H/2]$, as a function of the horizontal coordinate $x$ for $E=5\times 10^{-4}$, $St=1$ and $f=0$. The influence of the peripheral current is larger when the Rayleigh number is close to the critical Rayleigh number for flow in the bulk ((1.3)). Here, we have $Ra_{c}^{{bulk}}=6\times 10^{4}$, $Ra_{c}^{{wall}}=6.3\times 10^{4}$, giving (a) $Ra/Ra_{c}^{{bulk}}=1.66$, $Ra/Ra_{c}^{{wall}}=1.57$ and (b) $Ra/Ra_{c}^{{bulk}}=6.6$, $Ra/Ra_{c}^{{wall}}=6.3$.

3.1.3. Initial fluid layer height

The effective Rayleigh number at $t=0$ is determined by the initial height of the liquid $h_{0}$. In recent studies of convection-driven melting (e.g. Rabbanipour Esfahani et al. 2018; Favier et al. 2019), the initial liquid height is taken to be small fraction of the domain height $H$, such that $Ra_{{eff}}(t=0) < Ra_{c}$. Thus, convection begins only after an initial stage where melting occurs by the relatively slow diffusion of heat, which eventually leads to $Ra_{{eff}}>Ra_{c}.$ In our simulations with $h_{0}=H/2$, the initial $Ra$ is sufficiently large so that convection occurs immediately. While the melting history will obviously depend on $h_{0}$, this choice does not change the general conclusions drawn from our simulations. We show this in figure 11 by comparing the void area and number as a function of the height of the fluid layer in simulations with $h_{0}=1$ and $h_{0}=0.1$. Apart from initial transient differences, the curves follow very similar trajectories.

Figure 11. (a) The number of solid voids, and (b) the area of the solid voids as a function of the liquid height $h$, for $E=10^{-4}$, $Ra=10^{7}$, $Pr=5$, $St=1$, $f=0$.

3.1.4. Stefan number

Smaller Stefan numbers, as defined in (2.8), are associated with large latent heats and thus lead to lower melt rates (see e.g. Worster 2000), in which case simulations need to be run for longer times. However, the melting morphology we find is independent of Stefan number for the range studied ($St = 0.2\text {--}1$), which is shown by plotting the number and areas of the voids formed in figure 12. The same is found in the melting of pure solids driven by non-rotating convection (Rabbanipour Esfahani et al. 2018; Favier et al. 2019).

Figure 12. (a) The number of solid voids; and (b) the area of the solid voids as a function of the liquid height $h$, for $E=3.2\times 10^{-4}$, $Ra=2\times 10^{6}$, $Pr=5$, $f=0$.

3.2. Special cases

3.2.1. No-slip lower boundary: melting rates, flow structures and wall modes

In the simulations presented thus far, the fluid layer is bounded laterally and above by no-slip boundaries. Only the heated lower boundary is one of free slip. In rotating Rayleigh–Bénard convection, the role of the velocity boundary layers is as essential as in the non-rotating case (e.g. Rossby 1969; Liu & Ecke 2009; Schmitz & Tilgner 2010; Julien et al. 2012; King et al. 2012). Moreover, the critical Rayleigh number in (1.3) is largest for free-slip top and bottom boundaries, and smallest for one free-slip and one no-slip boundary; the case of two no-slip boundaries is intermediate between these cases (Chandrasekhar 1953; Boubnov & Golitsyn 1990). Despite this, for the parameter ranges considered here, the Nusselt number is larger for the case with no-slip upper and lower boundaries, owing to the interaction of the thermal and velocity boundary layers at the lower boundary (e.g. Rossby 1969). Thus, the melting rates are higher when the lower boundary is one of no slip as compared to one of free slip, as seen in figure 13.

Figure 13. The melting histories with either a no-slip or a free-slip lower boundary. The other parameters are identical, with $E=8\times 10^{-5}$, $Ra=7.8\times 10^6$, $Pr=5$, $St=1$, $f=0$. Due to the enhanced heat transport, the rate of melting is higher with a no-slip lower boundary.

Experiments show that columnar vortices in rotating convection show horizontally diffusive motion (see e.g. Noto et al. 2019). Because the phase boundary voids created by the heat transported through the columnar vortices are colocated, the latter can be arrested (and perhaps pinned) by the former. In our simulations, this effect is influenced by velocity boundary conditions, with horizontal motion suppressed in the case of no-slip boundaries. In figure 14, we show that the wall-modes that usually precess in a retrograde (i.e. clockwise as seen from above) direction are locked in place as the solid melts, an effect that is more prominent with a no-slip lower boundary than with a free-slip lower boundary.

Figure 14. A Hövmöller plot of the temperature, $\theta$, and the vertical velocity, $w$, at one of the vertical walls. The wall modes, which usually propagate clockwise, are locked in place once melting begins. The parameters are $E=8\times 10^{-5}$, $Ra=1.56\times 10^6$, $Pr=1$, $St=1$, $f=0$.

3.2.2. Horizontal periodicity

As we have seen, the presence of walls confining the flow in the horizontal directions leads to the generation of a peripheral current that can affect the melting of the solid. This peripheral flow is absent in a horizontally periodic system, as seen in figure 15. However, the columnar-vortical flow at $Pr=5$ and the resultant melt pattern reflecting the presence of these vortices, as well as the overall melt rate, both remain unchanged.

Figure 15. The solid–liquid interface (viewed from the solid side) at $t=400$ for $E=8\times 10^{-5}$, $Ra=7.8\times 10^{6}$, $Pr=5$, $St=1$, $f=0$, and (a) no-slip walls (b) periodic in the horizontal. The effects of the peripheral streaming flow seen in (a) as increased melting near the walls is absent in (b), although the voids and the overall rate of melting are very similar in the two cases.

3.2.3. Thermal diffusivity in the solid

The thermal diffusivity of the solid governs the amount of heat transported away from the solid–liquid interface and thus the melt rate (see (2.19)), with a larger diffusivity in the solid $\hat {\kappa }_{s}$ leading to smaller $u_m$. Figure 16(a) shows this effect for two values of $\hat {\kappa }_{s}=0.2$ and $\hat {\kappa }_{s}=5$, with $f=1$ (so that the upper boundary is at $\theta =-1$). For the largest value of the diffusivity, $\hat {\kappa }_{s}=5$, and the smaller melting rate (see, (2.12)), the horizontal drift of the columnar vortices is faster than the melt rate and hence we infer the vortices are not pinned in the voids. As a result, we see in figure 16(b) that the voids have smaller amplitudes.

Figure 16. (a) The melting histories for $E=8\times 10^{-5},$ $Ra=7.8\times 10^{6}$, $St=1$, $f=1$, $Pr=5$ with periodic boundary conditions in the horizontal for (i) $\hat {\kappa }_{s}=0.2$ and (ii) $\hat {\kappa }_{s}=5$. As $\hat {\kappa }_{s}$ increases the melt rate decreases. (b) The solid–liquid interface (viewed from the solid side) for $\hat {\kappa }_{s}=5$ at $t=500$. The system is periodic in the horizontal, and the other parameters are as in (a), but the voids are not as prominent.

3.3. Coupling of interfacial geometry and flow structure

We argued in § 3.2.3 that the phase boundary and the flow structures co-evolve, which is particularly well reflected in figure 3 showing the proportionality between the number and area of the vortices for $St=1$. Whilst we are unable to track individual vortices in our simulations, in figure 17 we assess their interaction with the voids by plotting the time evolution of the net vertical circulation, or vorticity, and the roughness, as characterized by the standard deviation of the liquid height $\sigma (h)$. We see that the rates at which both roughness and vorticity increase, decrease as the latent heat increases and that the roughness and the vorticity increase collinearly, which is a natural consequence of the conservation of potential vorticity. Indeed, we speculate that the increase in vorticity with latent heat shown in figure 17(b) is associated with the horizontal drift of the columnar vortices being faster than the evolution of the phase boundary. However, in order to assess such a scenario one must track individual vortices.

Figure 17. (a) The net circulation ($\varGamma = \iint \omega _z\,{\textrm {d}x}\,{\textrm {d}y}$) at $z=H/4$ and the roughness as a function of time; and (b) the net circulation $\varGamma$ as a function of the roughness, characterized by the standard deviation of the liquid height, $\sigma (h)$, for simulations with $E=3.2\times 10^{-4}$, $Ra=2\times 10^{6}$, $f=0$, $Pr=5$ and $St=0.05, 0.2, 1$ showing that the total vorticity and the roughness co-vary. The curves are computed using a running average over $10$ points, each spaced $20$ flow units apart.

3.4. Heat transport and the melting rate

In § 2.3 we noted that the initial and boundary conditions in most of the simulations reported here, except those in §3.2.3, are that the solid is at the melting temperature throughout, viz., $\theta (t=0)=0$, and the upper boundary is held at $\theta =0$. Therefore, the heat available for melting is transported by the fluid from the lower heated boundary to the solid and described by the integral form of energy conservation, (2.17), as

(3.2) \begin{align} &\rho\lambda (H/2)^2 U_{b}\left[\frac{\textrm{d}}{\textrm{d}t}\iiint(1-\chi)\,\textrm{d}x\,\textrm{d}y\,\textrm{d}z\right]\notag\\ &\quad = k_{l}\Delta T A^2 (H/2)\left[\left\langle -\frac{\partial\theta}{\partial z}\right\rangle _{z=0}\right]\nonumber\\ &\qquad - \rho C_{p}\Delta T (H/2)^{2}U_{b}\left[\frac{\textrm{d}}{\textrm{d}t}\iiint\theta\,\textrm{d}x\,\textrm{d}y\,\textrm{d}z\right], \end{align}

where the terms in square brackets are non-dimensional. Dividing by $k_{l}\Delta T A^2 H/2$ gives

(3.3)\begin{equation} {\frac{\left(Ra Pr\right)^{1/2}}{St}\frac{\textrm{d}h}{\textrm{d}t} =\left\langle -\frac{\partial\theta}{\partial z}\right\rangle _{z=0}- 2\left(Ra Pr\right)^{1/2}\frac{\textrm{d}\bar{\theta}}{\textrm{d}t},} \end{equation}

where

(3.4)\begin{equation} \bar{\theta}=\frac{1}{2A^2}\iiint\theta\,\textrm{d}x\,\textrm{d}y\,\textrm{d}z \end{equation}

is the average non-dimensional temperature over the simulation volume and

(3.5)\begin{equation} h=\frac{1}{A^2}\iiint(1-\chi)\,\textrm{d}x\,\textrm{d}y\,\textrm{d}z \end{equation}

is the volume-averaged dimensionless height of the fluid. The relative contributions of the sensible heating of the fluid and the melting of the solid to the heat balance are shown in figure 18. Initially, all the energy supplied to the system from the boundary heats up the liquid. For smaller $E$ and $Ra$, vertical motions are suppressed and hence so too is the delivery of the specific heat to the phase boundary, where melting may proceed (beginning here at about $t=50$). Once melting begins the latent heat draws down the sensible heat stored in the fluid and eventually a near steady balance between the energy delivered and that available for melting may be maintained. Hence, whilst the vigour of convection depends on $E$ and $Ra$, such a balance between the heat input at the lower boundary and the latent heat of fusion requires quasi-steady rotating convection.

Figure 18. The terms in (3.3), with $Pr=5$ and $f=0$ for (a,c) $E=10^{-4}$ $Ra=10^{7}$; (b,d) $E=3.2\times 10^{-4}$, $Ra=2\times 10^{6}$. The Stefan numbers are (a,b) $St=1$; (c) $St=0.1$; (d) $St=0.05$. Note that the quasi-steady state of convection in the fluid described by (3.3) breaks down when the voids in the solid reach the upper boundary and fluid comes into direct contact with the container surface at $t=340$ in (b).

We see in figure 18(b) that the quasi-steady state of convection in the fluid described by (3.3) breaks down at $t=340$ when fluid comes into contact with the upper solid boundary through the voids in the solid. Note that the slight mismatch between $\left \langle -({\partial \theta }/{\partial z})\right \rangle _{z=0}$ and the sum

(3.6)\begin{equation} \frac{(Ra Pr)^{1/2}}{St}\frac{\textrm{d}h}{\textrm{d}t}+2(Ra Pr)^{1/2}\frac{\textrm{d}\bar{\theta}}{\textrm{d}t} \end{equation}

in figure 18 is a consequence of the coarse time discretization used in calculating the time derivatives in the plots.

Additionally figure 18 shows that, when the specific heat stored in the convecting fluid is small, i.e. when the Stefan number is small, there is a nearly steady balance between the heat supplied at the base of the cell and the melt rate. As the fluid interior cools slightly in time this is balanced by a slight increase of the melt rate and the heat input from the lower boundary, as seen in figure 18(c,d). The temperature in the liquid is, of course, not uniform in space. Indeed, as shown in figure 19, the structure of the mean temperature gradient in the fluid is reminiscent of non-rotating high $Ra$ convection, with a thermal boundary layer at the base and a nearly isothermal interior. However, the phase change at the ramified upper boundary maintains the average temperature near the melting point. This situation can be treated by approximating (3.3) using only the first two terms, viz.,

(3.7)\begin{equation} \frac{\left(Ra Pr\right)^{1/2}}{St}\frac{\textrm{d}h}{\textrm{d}t} =\frac{Nu}{h}, \end{equation}

where $Nu$ is the Nusselt number – the total heat flux scaled by the conductive heat flux – across the fluid region.

Figure 19. The area-averaged solid fraction $\langle \chi \rangle$ and temperature $\langle \theta \rangle$ as a function of the vertical coordinate $z$ at $t=400$, for the case $E=10^{-4}$, $Ra=10^{7}$, $Pr=5$, $St=1$ and $f=0$.

In § 3.1 we showed that for most combinations of parameters examined here, the phase boundary is ramified, so that the solid depth varies substantially in the horizontal. In consequence, we see from figure 19 that within the broad average transition region from fluid to solid the average temperature relaxes to the bulk melting temperature. Therefore, we take the domain-averaged $h$ (see also § III of Rabbanipour Esfahani et al. 2018) when considering the quasi-steady balance in (3.7). We note, however, that we understand that there are three-dimensional heat fluxes in the interfacial region, which are simpler to treat when the phase boundary has small amplitude variations, such as in the non-rotating case (e.g. Favier et al. 2019; Toppaladoddi & Wettlaufer 2019). Another perspective is that for a vortex-induced highly ramified interface, the interfacial region might be considered as a ‘mushy layer’, as observed in binary systems (Worster 2000), wherein there is two-phase, two-component coexistence and the condition of marginal equilibrium holds. Clearly here there are no impurities, but we can see in figure 19 the relaxation towards equilibrium of the average temperature and enthalpy through the mixed phase region.

For geostrophic convection, the average $Nu$ can be expressed in terms of the Rayleigh number and the critical Rayleigh number, using (2.20), (2.21) and (1.3), as

(3.8)\begin{equation} Nu=C\left(\frac{Ra}{Ra_{c}^{{bulk}}}\right)_{{eff}}^{\beta}, \end{equation}

where $\beta$ is in general a function of $(Ra/Ra_{c}^{{bulk}})_{{eff}}$ and $C$ is a numerical prefactor that may depend on $Pr$. For large values of $Ra/Ra_{c}^{{bulk}}$, two values have been suggested in the literature; $\beta =3$ (Boubnov & Golitsyn 1990; King et al. 2012) and $\beta =3/2$ (Julien et al. 2012), the latter finding $C=(1/25)Pr^{-1/2}$. For more modest values of $Ra/Ra_{c}^{{bulk}}$, Ravichandran & Wettlaufer (2020) found $\beta =3/4$ and Liu & Ecke (2009) found $\beta =2/7$. In the limit of large $Ro$, that is in the classical non-rotating Rayleigh–Bénard convection regime, one finds, with a different prefactor than in (3.8), $\beta =1/3$ up to $Ra=10^{15}$ (Doering, Toppaladoddi & Wettlaufer 2019; Doering 2020a,b; Iyer et al. 2020).

In figure 20 we plot the Nusselt number, calculated using (3.7), vs the effective Rayleigh number as melting proceeds. In the quasi-steady state the curves for different $St$ collapse with $E$ and $Ra$ dependent slopes, suggesting that although (3.8) provides an ideal organizing principle for our simulations, we are unable to determine the associated exponent given our parameter range (see e.g. Stumpf & Porter 2012).

Figure 20. The Nusselt number plotted as a function of $(Ra/Ra_{c})_{{eff}}$, calculated using the mean fluid height $h(t)$ and (2.20) and (2.21).

3.5. Maximal phase boundary roughness and maximal heat flux

We conclude § 3 with the observation that the roughness of the phase boundary continuously increases and reaches a maximum approximately simultaneously with the Nusselt number. As seen in figure 18, the heat supplied at the bottom boundary, and the melt rate of the solid, are approximately independent of time and hence the left-hand side of (3.7) is approximately constant. Therefore, the Nusselt number increases linearly with the liquid height $h$ and reaches a maximum when the voids in the solid reach the upper boundary of the cell. We again characterize the roughness using the standard deviation of the liquid height, $\sigma (h)$, which we observe reaches a maximal value when the voids reach the upper boundary, namely when there is fluid in contact with the upper boundary. Further melting reduces the roughness. The correlation between $Nu$ and $\sigma (h)$ is shown in figure 21(a), where we see that the maximal Nusselt numbers are reached before the roughness of the solid–liquid interface becomes maximal, with the interval between the maxima increasing as the Stefan number decreases (and the melt rate decreases). Smaller Stefan numbers lead to voids of unequal depths, with some voids reaching the upper boundary before others. The decrease of the interface roughness associated with the former is compensated, for a limited period, by the continued deepening of the latter.

Figure 21. The Nusselt number $Nu$ and standard deviation of the phase boundary height, $\sigma (h)$, plotted as a function of time for $f=0$ and (a) $E=3.2\times 10^{-4}$, $Ra=2\times 10^6$, $Pr=5$, $h_{0}=1.0$, for a range of $St$; and for (b) $St=1$ with two combinations of $E,Ra$ and $h_{0}=0.1$. The correlation between $Nu$ and $\sigma (h)$ is evident in all of these cases. In (c), we show the roughness data in (a) with the time coordinate rescaled by the Stefan number. Thus, we see the Nusselt number maxima occurring at smaller $t \times St$ for smaller $St$, from which we expect that the data for $St=0.05$ will follow this trend and reach a maximum, were we able to run longer simulations in that case.

Since the areas and number of voids depend on the flow parameters (figure 5), the maximum value of $\sigma (h)$ depends on these parameters as well, with the thinner vortices in flows with smaller $Ro_{c}$ ((2.22)) leading to narrower voids and thus a rougher interface (see e.g. figure 21b). In particular, the continued increase of $\sigma (h)$ in figure 21(b), where the initial liquid height $h_{0}=0.1$, shows that the voids formed by the columnar vortices will continue to penetrate deeper into the solid with time, only being limited by the depth of the solid itself. The curves for $St=0.05$ in figure 21(a) have not reached their maxima. However, the rescaling of the data in figure 21(c) suggests that the time interval between the maximal $Nu$ and the maximal $\sigma (h)$ will further increase for $St=0.05$, and we expect that a maximum will be reached were we able to run longer simulations in that case.

In non-rotating turbulent Rayleigh–Bénard convection, with Dirichlet boundary conditions and periodically rough boundaries, Toppaladoddi, Succi & Wettlaufer (2015, 2017) showed that, for a given roughness wavelength, there is a ratio of the thermal boundary layer thickness to the roughness amplitude that optimizes $Nu$. Moreover, this enhancement of heat transport is a general consequence of roughness, observed for a wide range of geometries from rough on all surfaces to fractal boundaries (Roche et al. 2001; Goluskin & Doering 2016; Toppaladoddi et al. 2020). In these situations, the systems are in a statistical steady state. Here, with the geometry free to evolve subject only to the underlying conservation laws, both the roughness of the phase boundary and the Nusselt number increase with time as the solid phase melts according to (3.7). The Stefan number dependence of the observed correlation between the vorticity and the interfacial roughness shown in figure 17 underlies this process.

4. Conclusion

We have studied the melting of a pure solid by the convection of its liquid phase when the former overlies the latter and the entire system rotates about an axis parallel to gravity. The width of the system is twice its depth and we have examined ranges of the Ekman, Rayleigh and Prandtl numbers predominantly corresponding to moderately rotating Rayleigh–Bénard convection.

There are three regimes of flow that influence the morphology of the phase boundary. First, when the Rayleigh number is greater than the bulk critical value, $Ra>Ra_{c}^{{bulk}}$ ((1.3)), the flow takes the form of columnar vortices. Second, in confined geometries there is a streaming flow close to the lateral walls of the container. This occurs when $Ra_{c}^{{wall}} < Ra < Ra_{c}^{{bulk}}$, where $Ra_{c}^{{wall}}$ is given by (1.4) (Herrmann & Busse 1993). Third, in the periodic geometry, there is no flow for $Ra < Ra_{c}^{{bulk}}$. We found that the number of melt voids in the solid is proportional to the number of heat transporting vortices present, which in turn increases as the convective Rossby number decreases and rotational effects become dominant. We showed that the overall melting rate is a non-trivial function of the flow parameters; for the same $Ra/Ra_c$, melting rates are smaller for larger Prandtl numbers and smaller $E$. Moreover, we found that the phase boundary morphology can be highly ramified or relatively smooth, reflecting the nature and number of rotationally controlled vortices transporting heat across the evolving fluid layer. Lastly, we showed that the peripheral streaming current characteristic of rotating Rayleigh–Bénard convection may become ‘locked’ in place due to the coupling between the flow and the melting of the solid.

For large values of the latent heat of fusion, characterized by the Stefan number, we found a quasi-steady geostrophic convective state in which the net vertical heat flux is nearly constant over long time intervals. This leads to a situation in which the constant heat supplied at the base balances the melt rate. In the case of non-rotating binary systems, it is now well known that the fluid mechanics of solidification lead to complex phase boundary geometries and their associated transport phenomena (e.g. Huppert 1990; Sullivan, Liu & Ecke 1996; Worster 2000; Philippi et al. 2019). Here, in contrast, in a pure system, we find that convective and rotationally controlled vortices alone can create ramified phase boundaries. While no obvious optimization of the Nusselt number is seen as a consequence of the increasing boundary roughness, that roughness evolves in time in a unique manner coupled to the rotationally influenced evolving buoyancy of the liquid phase. The associated void structure in the solid will affect the mechanical and thermal properties of materials formed in such circumstances. Thus, the inclusion of compositional effects with the rotational processes studied here will open a new set of questions regarding the structure of partially molten rotating systems. Finally, we note that in astrophysical and geophysical problems wherein rotational effects are important, assumptions of planarity of the phase boundary should therefore be made with care.