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 Reference Davis2001). 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 Reference Huppert1990; Worster Reference Worster2000). In astrophysics, they underlie planet formation (e.g. Armitage Reference Armitage2020), wherein for example the proto-Earth was believed to rotate approximately ten times faster than today (e.g. Cuk & Stewart Reference Cuk and Stewart2012), and the growth of neutron star crusts (e.g. Baym et al. Reference Baym, Hatsuda, Kojo, Powell, Song and Takatsuka2018), 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
\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 Reference Chandrasekhar1961).
In direct analogy with stratification in non-rotating systems, rotation suppresses vertical motions due to buoyancy (Veronis Reference Veronis1970). Therefore, the critical Rayleigh number above which convection occurs is a function of the rotation rate of the system (Chandrasekhar Reference Chandrasekhar1953). The Ekman number is the relevant non-dimensional rotation rate and is
\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 Reference Chandrasekhar1953), 
$Ra_{c}$ is
\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 Reference Zhong, Ecke and Steinberg1991; Ecke, Zhong & Knobloch Reference Ecke, Zhong and Knobloch1992) is, in the limit of large 
$E^{-1}$, given by (Herrmann & Busse Reference Herrmann and Busse1993)
\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 Reference Horn and Schmid2017; De Wit et al. Reference De Wit, Guzmán, Madonia, Cheng, Clercx and Kunnen2020; Favier & Knobloch Reference Favier and Knobloch2020; Zhang et al. Reference Zhang, Van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020), even when there are severe obstacles in the way (Favier & Knobloch Reference Favier and Knobloch2020). 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. Reference De Wit, Guzmán, Madonia, Cheng, Clercx and Kunnen2020).
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 Reference Boubnov and Golitsyn1986, Reference Boubnov and Golitsyn1990; Zhong et al. Reference Zhong, Ecke and Steinberg1991; King et al. Reference King, Stellmach, Noir, Hansen and Aurnou2009; Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015). These vortices are predominantly cyclonic near the upper and lower boundaries, with equal numbers of cyclonic and anticyclonic vortices in the interior (Boubnov & Golitsyn Reference Boubnov and Golitsyn1986; Zhong et al. Reference Zhong, Ecke and Steinberg1991; Vorobieff & Ecke Reference Vorobieff and Ecke1998; Kunnen, Clercx & Geurts Reference Kunnen, Clercx and Geurts2010), thereby transporting heat from the boundaries (Sakai Reference Sakai1997). 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 Reference Boubnov and Golitsyn1990; King, Stellmach & Aurnou Reference King, Stellmach and Aurnou2012; Shi et al. Reference Shi, Lu, Ding and Zhong2020), 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
\begin{gather} \rho_{s} =\rho_{l}\left(=\rho\right), \end{gather}
\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
\begin{gather} \mathcal{H}_{s} =\rho C_{p}T,\text{ and } \end{gather}
\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
\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
\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
\begin{equation} \theta = \frac{T-T_{m}}{\Delta T} \end{equation}is defined to be the non-dimensional temperature, and
\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
\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
\begin{equation} \chi=1-\max\left[0,\min\left(1,St\, \phi\right)\right], \end{equation}and hence the temperature follows as
\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. Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018; Favier, Purseed & Duchemin Reference Favier, Purseed and Duchemin2019), 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
\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
\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}
\begin{gather}\frac{\textrm{D}\theta}{\textrm{D}t} = \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\kappa\boldsymbol{\nabla}\theta\right), \text{ and } \end{gather}
\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
\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}
\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}
\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
\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
\begin{gather} Ra_{{eff}} = Ra\left[\frac{h\left(t\right)}{H/2}\right]^{3},\quad\text{and} \end{gather}
\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
\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 Reference Cushman-Roisin and Beckers2011, 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 (Reference Ravichandran and Wettlaufer2020) 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. Reference Diwan, Prasanth, Sreenivas, Deshpande and Narasimha2014; Prasanth Reference Prasanth2014; Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2020; Ravichandran, Meiburg & Govindarajan Reference Ravichandran, Meiburg and Govindarajan2020). 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. (Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018) 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 (Reference Kevlahan and Ghidaglia2001), 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
\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. Reference Zhong, Ecke and Steinberg1991; Sakai Reference Sakai1997; Vorobieff & Ecke Reference Vorobieff and Ecke1998), 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 Reference Popinet2018). 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 Reference Stumpf and Porter2012) 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. Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018; Favier et al. Reference Favier, Purseed and Duchemin2019), 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 Reference Worster2000), 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. Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018; Favier et al. Reference Favier, Purseed and Duchemin2019).
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 Reference Rossby1969; Liu & Ecke Reference Liu and Ecke2009; Schmitz & Tilgner Reference Schmitz and Tilgner2010; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012; King et al. Reference King, Stellmach and Aurnou2012). 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 Reference Chandrasekhar1953; Boubnov & Golitsyn Reference Boubnov and Golitsyn1990). 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 Reference Rossby1969). 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. Reference Noto, Tasaka, Yanagisawa and Murai2019). 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
\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
\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
\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
\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
\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.,
\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. Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018) 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. Reference Favier, Purseed and Duchemin2019; Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019). 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 Reference Worster2000), 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
\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 Reference Boubnov and Golitsyn1990; King et al. Reference King, Stellmach and Aurnou2012) and 
$\beta =3/2$ (Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012), the latter finding 
$C=(1/25)Pr^{-1/2}$. For more modest values of 
$Ra/Ra_{c}^{{bulk}}$, Ravichandran & Wettlaufer (Reference Ravichandran and Wettlaufer2020) found 
$\beta =3/4$ and Liu & Ecke (Reference Liu and Ecke2009) 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 Reference Doering, Toppaladoddi and Wettlaufer2019; Doering Reference Doering2020a,bReference Doering; Iyer et al. Reference Iyer, Scheel, Schumacher and Sreenivasan2020).
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 Reference Stumpf and Porter2012).
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 (Reference Toppaladoddi, Succi and Wettlaufer2015, Reference Toppaladoddi, Succi and Wettlaufer2017) 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. Reference Roche, Castaing, Chabaud and Hébral2001; Goluskin & Doering Reference Goluskin and Doering2016; Toppaladoddi et al. Reference Toppaladoddi, Wells, Doering and Wettlaufer2020). 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 Reference Herrmann and Busse1993). 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 Reference Huppert1990; Sullivan, Liu & Ecke Reference Sullivan, Liu and Ecke1996; Worster Reference Worster2000; Philippi et al. Reference Philippi, Berhanu, Derr and du Pont2019). 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.
Funding
Computational resources from the Swedish National Infrastructure for Computing (SNIC) under grants SNIC/2018-3-580, SNIC/2019-3-386 and SNIC/2020-5-471 are gratefully acknowledged. Computations were performed on Tetralith. The Swedish Research Council, under grant no. 638-2013-9243, is gratefully acknowledged for support.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation of the enthalpy method
We validate the enthalpy method used here by comparing the numerical solution to the one-dimensional analytical solution for a purely conducting case (e.g. Worster Reference Worster2000). We then study the convergence of the method with grid resolution in a case with fluid convection.
A.1. Melting by conductive heat transfer
Consider a semi-infinite solid layer in the region 
$z>0$ at the melting temperature. The boundary at 
$z=0$ is held at 
$\theta =1$. The solid melts, forming a liquid layer of height 
$h(t)$ given by
\begin{equation} h=2\xi\sqrt{\kappa t}, \end{equation}
where 
$\xi$ is the solution of the transcendental equation deriving from the Stefan condition,
\begin{equation} \xi\exp(\xi^{2})\text{erf}(\xi)=\frac{St}{\sqrt{\rm \pi}}. \end{equation}
In figure 22(a) the analytical solution of the Stefan problem is compared with a numerical solution of (2.17) in one dimension with the boundaries at 
$z=0$ and 
$z=H=0.5$. Next, we consider a case where there is already some liquid (at 
$\theta =0$) present in the region 
$0 < z < z_{0}=0.05$, with solid at the melting temperature in the region 
$z_{0} < z < H$ at 
$\theta =0$. The boundaries are held at 
$\theta (z=0)=1$ and 
$\theta (z=H)=0$. The numerical solution in one dimension is compared with the solution from the three-dimensional solver, and the amount of unmelted solid plotted as a function of time in figure 22(b). In both these cases, the one-dimensional solution is obtained using fourth-order Runge–Kutta integration; the three-dimensional solver uses a second-order Adams–Bashforth scheme (as described in § 2.4).
Figure 22. (a) The liquid height 
$h(t)$ from the one-dimensional numerical solution and the analytical solution of the Stefan problem ((A1)). In the numerical solution 
$h(t)$ is bounded by the height of the domain, 
$0.5$. (b) For the alternate initial conditions (see text), the amount of unmelted solid from the numerical solution from the finite-volume solver is compared with the analytical solution in one dimension. The parameters are 
$\kappa =0.01$, 
$St=1$.
For the single-component, two-phase systems considered here, the solid–liquid interface is sharp. In the numerical simulations, this interface is defined as the region where 
$0<\chi <1$, and is distributed over a finite number of grid points. This is shown in figure 23 where the mask 
$\chi$ and the temperature 
$\theta$ are plotted on a vertical line through the peak of the void in the solid region. The mask function 
$\chi$ varies from 
$0$ to 
$1$ over a distance of about 
$\delta z=0.008$, which is 
$2$ grid points in the 
$256^{2}\times 128$ simulations. This is similar to results obtained by Couston et al. (Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2021), and those prescribed (in their formulation) by Favier et al. (Reference Favier, Purseed and Duchemin2019) who use a nominal interface thickness of half the grid spacing. We note that for the range of values of 
$\eta$ used here, the thinness of the interface is not affected by changes in the grid resolution or in the penalization parameter, as seen from figure 23(b), with 
$\eta =10^{-3}$.
Figure 23. The variation of the solid mask and the temperature through the solid–liquid interface. Parameters: 
$E=8\times 10^{-5},$ 
$Ra=7.8\times 10^{6}$, 
$Pr=5$, 
$St=1$, 
$f=0$, grid spacing 
$\textrm {d}z\approx 0.004$ with (a) 
$\eta =2\times 10^{-3}$, (b) 
$\eta =10^{-3}$. The dotted lines in (a) show the average thickness of the thermal boundary layer at the heated lower boundary.
A.2. Melting by convective heat transfer
The grid dependence of the accuracy of our solution method is examined as follows. We use the geometry in Appendix A2 of Favier et al. (Reference Favier, Purseed and Duchemin2019), and 
$Ra=1.25\times 10^5$, 
$Pr=1$ and 
$St=1$, with an initial temperature perturbation of
\begin{equation} \theta(t=0) = 1 - z + A \sin(2 {\rm \pi}x) \sin({\rm \pi} z). \end{equation}
The resulting velocity and temperature fields at 
$t=56$ are plotted in figure 24. The location of the solid–liquid interface is given by the liquid height 
$h$ from (3.5), and is plotted as the grid resolution is varied in figure 25(a). We then use the solution at the highest grid resolution (
$N=1024$) as a reference, and linear interpolation to find the interface location at intermediate points. The root-mean-square error is plotted as a function of 
$N$ in figure 25(b), showing that the error decreases as 
$N$ is increases, with an exponent between 
$1$ and 
$2$, as also reported by Favier et al. (Reference Favier, Purseed and Duchemin2019).
Figure 24. (a) The vertical velocity and (b) the temperature fields at 
$t=56$ for simulations at the highest resolution (
$N=1024$), with 
$Ra=1.25\times 10^5$, 
$Pr=1$ and 
$St=1$.
Figure 25. (a) The liquid height as the grid resolution is varied, and (b) the root-mean-square error as a function of the resolution.
Appendix B. Penalization parameter
The volume penalization method has a tuneable parameter 
$\eta$. The principle of the volume penalization method is to treat the solid as a porous medium of vanishing porosity. The use of a finite value for 
$\eta$ creates a velocity boundary layer of size 
$(\nu \eta )^{1/2}$ in the solid. Engels et al. (Reference Engels, Kolomenskiy, Schneider and Sesterhenn2015) showed that the optimal value of 
$\eta$ is such that the grid spacing is comparable to the boundary layer thickness, namely 
${\textrm {d}x}\sim (\nu \eta )^{1/2}$. All of our results are reported with the penalization parameter 
$\eta =2\times 10^{-3}$ (§ 2.4), satisfying this requirement.
In detail the melting process is influenced by the boundary layer and hence depends on 
$\eta$. As seen in figure 26, upon reduction of 
$\eta$ by a factor of 
$2$, the melt rate changes by only a few per cent. Therefore, the latent heat flux and the quasi-steady balance described by (3.7) underlying the results shown in figures 18(b) and 21 are insensitive to the choice of 
$\eta$. Snapshots of the interface shown in figure 27, and the plots of the number and areas of the voids shown in figure 28 demonstrate the persistence of the central behaviour; convective vortices etch voids into the solid, and the number of voids are proportional to the number of vortices. Thus, as noted in § 3.5, 
$Nu(t)$ and the maximal interface roughness depend on 
$\eta$, but the correlation between 
$Nu$ and 
$\sigma (h)$ shown in figure 21 do not.
Figure 26. (a) The melting history and (b) the melting Nusselt number for 
$E=8\times 10^{-5},$ 
$Ra=7.8\times 10^{6}$, 
$Pr=5$, 
$St=1$, 
$f=0$, as 
$\eta$ is varied. The difference in the total amount of solid melted changes by only approximately 
$5\,\%\text {--}10\,\%$ over 
$250$ flow units when 
$\eta$ is halved from 
$2\times 10^{-3}$ to 
$10^{-3}$. The differences in the melting rates are even smaller. As a result, the Nusselt number also changes by only approximately 
$5\,\%\text {--}10\,\%$ as the 
$\eta$ is halved.
Figure 27. Snapshots of the phase boundary at 
$t=500$ for the case 
$E=8\times 10^{-5},$ 
$Ra=7.8\times 10^{6}$, 
$Pr=5$, 
$St=1$, 
$f=0$, with (a) 
$\eta =2\times 10^{-3}$ and (b) 
$\eta = 10^{-3}$. The number and area of the voids, as well as the overall amount of melting (noting that the figures are plotted at the same time 
$t=500$), can be seen to be insensitive to the penalization parameter.
Figure 28. The void areas and the number of voids formed in the solid for different values of the penalization parameter. The number and area of the voids can be seen to be insensitive to the penalization parameter.
Open access
