Melting driven by rotating Rayleigh-B\'enard convection

We study the melting of a horizontal layer of a pure solid above a convecting layer of its fluid rotating about the vertical axis using numerical methods. In the rapidly rotating regime, and for the Rayleigh numbers of order $10^7$ considered here, convection takes the form of columnar vortices. Since these vortices transport heat from the bottom surface to the upper boundary, the melt pattern reflects the number and size of the columnar vortices, which in turn depend on the Prandtl, Reynolds, Rossby and Stefan numbers of the system, and on whether we treat periodic or confined horizontal geometries. The phase boundary can be highly ramified, reflecting the nature and number of heat transporting vortices. Whereas the number of vortices and the melt regions they produce increase with Reynolds number, the average area of each vortex decreases and hence so too does the average melt rate. In addition to the Stefan number, the overall melt rate also depends on the velocity boundary condition on the lower boundary. For large values of the latent heat of fusion, a quasi-steady geostrophic convective state is reached in which the net vertical heat flux, or Nusselt number, reaches nearly constant maximal values over long time intervals, so that the constant heat supplied at the base balances the melt rate. Commensurate with this, we find that the interfacial roughness is also maximal, independent of the flow parameters. The confluence of processes responsible for the range of phase boundary geometries found should influence the treatment of moving boundary problems in mathematical models, particularly those in astrophysical and geophysical problems where rotational effects are important.


I. INTRODUCTION
The coupling between a solid and the liquid from which is 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 [e.g., 1]. Whereas, for example, when a pure melt is solidified from above in the absence of convection the phase boundary is planar, convective flow patterns lead to non-planar interfaces. However, the uncontrolled interplay of convection, rotation, and phase change determines the dynamics of many geophysical and astrophysical systems. Indeed, such processes operate from Earth's core to the principal components of the cryosphere [e.g., 2,3]. In astrophysics, they underlie planet formation [e.g., 4], wherein for example the proto-Earth was believed to rotate about ten times faster than today [e.g., 5], and the growth of neutron star crusts [e.g., 6], amongst many other phenomena. The confluence of dynamic and thermodynamic processes in such systems is highly complex and involves multiple timescales, components and phases.
Here, we study a simplified system of a singlecomponent rotating phase boundary heated from below. Initially the lower half of the container is liquid and the upper half is solid. Rotationally controlled thermal buoyancy brings heat to the upper boundary to drive phase change. The morphology of the melting solid results from the nonlinear interactions between convection and melting at the solid-liquid interface.
The paper is organized as follows. We describe the structure of the problem in §II, providing details of the phase change treatment used; the approximations made; the relevant physical scales and the nondimensionliza-tion; the boundary and initial conditions; and the numerical algorithm used to solve the governing equations. In §III, we discuss the effects of the control parameters of the problem on the phase boundary morphology and in §IV, we analyze the results to obtain the melting rates and the Nusselt number as a function of the problem parameters.

II. STRUCTURE OF THE PROBLEM
Our study geometry is a cube of dimensions L×L×H, with gravity g in the −z direction, and rotating about the +z axis with an angular velocity Ω. The aspect ratio is taken to be A = L/H = 2. The upper half of the domain (z/H ≥ 1/2) is initially solid, and the system is heated from below with a constant temperature difference between the upper and lower boundaries. The heat transported by thermal convection melts the solid in the upper half of the domain, and we study how this melting proceeds with time.

A. Enthalpy Method
We employ a mixture theory approach to tracking the solid region, such that a solid fraction variable χ varies from 0 in the liquid state to 1 in the solid state. We describe the system in terms of its enthalpy, following [7] and the derivation in [8]. The heat capacities of the solid and liquid phases are C s and C l respectively, and the latent heat of fusion is L f . The solid and liquid enthalpies are H s = ρ s C s T, and (1) respectively, where ρ s and ρ l are the associated densities. The enthalpy of the solid phase at the melting temperature T m is H 0 = ρ s C s T m , and that of a mixture of solid and liquid phases with solid volume fraction χ is given by We nondimensionalize the enthalpy as where ∆T is the difference between the temperature of the lower boundary and the melting temperature. Thus, with which is a modified latent heat of fusion. Thus, where is the nondimensional temperature, and is the Stefan number. Note that the Stefan number is also sometimes defined as the inverse of the definition used here. We note that in the purely solid state, χ = 1 and θ ≤ 1, so that φ ≤ 0. In the purely liquid state, χ = 0 and φ = θ + St −1 . In the mixed phase, 0 < χ < 1 and θ ≡ 0 by definition. Thus, the equation of state (8) can be inverted to give the solid fraction in terms of the enthalpy as and the temperature as Here, for simplicity, we only consider the case where with ρ and C p being constants, so that Thus, in a pure solid, θ = φ, and in a pure liquid, θ = φ − St −1 . In a system with a single pure substance χ = 0 or 1 everywhere, except in the vicinity of the phase boundary, which needs to be thin [see e.g., 9,10] and our simulations obey this requirement. 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 where (∇T ) s and (∇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.

B. Governing Equations
The equations of motion that govern the evolution of the velocity u, and the enthalpy φ, defined in Eq. (8), are as follows. We study the rotating Oberbeck-Boussinesq equations with the assumptions in Eqs. (13) and (14), which are where α is the coefficient of thermal expansion, ν is the kinematic viscosity of the fluid, and κ = χκ s + (1 − χ) κ l is the local thermal diffusivity. These equations are nondimensionalised using the temperature scale ∆T of Eq. (8), and the horizontal extent of the domain, L, as the length scale, giving a buoyancy velocity as U b = (gα∆T L) 1/2 . Using these scales, the dimensionless equations of motion become where Re = U b L/ν is the Reynolds number, P r = ν/κ l is the Prandtl number, Ro = U b /(ΩL) is the Rossby number, andκ = κ/κ l is the ratio of the local thermal diffusivity to the diffusivity in the liquid. The Stefan condition (Eq. 16) in nondimensionalized form is given by whereκ s = κ s /κ l is the nondimensional thermal diffusivity in the solid. The Rayleigh and Ekman numbers are respectively. Finally, in the solid there is only heat conduction and hence u = 0 in Eqs. (20)(21).

C. Initial and Boundary conditions
At t = 0, the bottom half of the domain is liquid, with enthalpy φ = St −1 ; the top half has enthalpy φ = 0. The upper and lower boundaries are held at temperatures θ = −f and θ = 1 respectively (f = 0 except in §III F). The lateral boundaries are insulating, no-slip walls. No-slip conditions are also applied at the freely evolving phase boundary, at which the temperature is θ = 0. A free-slip velocity condition applies on the lower boundary, except in §III D where we study the influence of the no-slip velocity boundary condition.

D. Numerical Simulations
Equations (20 -21), together with Eq. (11), are solved using the finite volume solver Megha-5 on a uniform grid in all three space directions [11][12][13][14]. The requisite velocity conditions in the resulting arbitrarily shaped solid region are applied using the volume-penalization method of [15], wherein the solid is modeled as a porous medium with vanishing porosity. This amounts to adding a term − χ η u to the right hand side of Eq. (20), where η 1 is the penalization parameter.
The numerical implementation of the enthalpy method used here is validated through comparison with the analytical solution of the Stefan problem for the purely conducting case, as shown in Appendix A. Our simulations are performed with grids of 512 2 × 256 gridpoints for Re = 10 4 and 256 2 × 128 for Re = 5000, with a timestep of δt = 5 × 10 −4 . The results presented are independent of the grid resolution. We use a penalization parameter of η = 10 −3 , and discuss the effects of changing η in Appendix B.

III. MELTING MORPHOLOGY
In the flow studied here convection is influenced by rapid rotation. In such flows, the critical Rayleigh number below which there is no flow in the bulk, for one free-slip one no-slip boundary each (and periodic boundary conditions in the horizontal), in the limit of large E −1 [16], is Below this Rayleigh number, the balance between Coriolis and pressure gradient effects, or the geostrophic balance, suppresses vertical gradients in the flow, and, for the initial conditions employed here, the fluid remains in solid-body rotation and the velocity is zero in the rotating frame.
In confined geometries, however, finite flow can be seen for Rayleigh numbers lower than Ra bulk c , but higher than a critical Rayleigh number for the so-called 'wall-mode' [17,18] which, in the limit of large E −1 , is [19] Ra wall Thus, for the majority of cases we report here, there is no flow for Ra < Ra wall c ; the flow takes the form of a streaming flow close to the lateral walls of the container for Ra wall c < Ra < Ra bulk c ; columnar vortices are seen for Ra > Ra bulk c . The nature of the flow, as we show below, determines the morphology of the phase boundary.
We note that for the single component, two-phase system considered here, this interface has to be sharp; i.e. regions with 0 < χ < 1 are not strictly permitted. Numerically, as we note in §II A, χ varies smoothly from 0 to 1 over a finite number of gridpoints (see Fig. 18 in Appendix A). For the purposes of plotting, the solid-liquid interface is taken to be the iso-surface χ = 0.5, although as should be clear from above, any value between 0 and 1 would give the same results.
A. Influence of the Prandtl number Figure 1 shows the morphology of the solid as it melts. For P r = 5, the flow takes the form of columnar vortices spanning the entire liquid height. On the other hand, for P r = 1 a peripheral streaming flow of alternating upand downwelling regions results. In both cases, the flow structures impinge on the solid from below, and melt the solid in patterns that reflect these structures. Figure 1 shows the resultant solid-liquid interface at t = 150. The solid liquid interface, which is at z = H/2 = 0.25 at t = 0, shifts upwards to z > 0.25. When columnar vortices are present, the solid melts in an array of inverted lobes, the structure of which reflect the vortex properties. Also shown in Fig. 1 are vertical cross-sections of the temperature distribution on the midplane defined as y = 0. These show the influence of the peripheral flow in the enhanced melting that occurs near the lateral walls. These effects are also seen in Fig. 2. This can be quantified by plotting the azimuthally averaged height of the solid-liquid interface, as shown in Fig. 7.

B. Influence of the Rossby number
As the Rossby number is increased, rotation becomes less dominant, and the P r-dependent classical convective vortex structures [20] become less prominent. The increase in Ro in Fig. 2 relative to Fig. 1, all other parameters being the same, has two consequences. First, for P r = 5, (Figs. 1a and 2a), the number of columnar vortices decreases (see also Fig. 10a) and the vortices become wider (Fig. 8b). Second, increasing Ro increases the heat transfer and the overall rate at which the solid melts (compare Figs. 14(a,c) and 14(b,d)). Furthermore, for P r = 1, and Ro = 0.1 in Fig. 1   At this higher Ro, the columnar vortices seen for P r = 5 are wider and fewer in number than in Fig. 1, as reflected in the morphology of the solid-liquid interface. For P r = 1, the influence of the peripheral streaming is still evident, but, owing to the nature of the flow, the interface is no longer flat (see also §IV A). These should be compared with the corresponding cases in Fig. 1.

C. Influence of the Reynolds number
The Reynolds number (Re = U b L/ν) quantifies the balance between the thermal forcing and the viscosity. Fig. 3 shows the solid-liquid interface for Ro = 0.2, Re = 5 × 10 3 , and may be compared with Fig. 1. For P r = 1, the four-fold decrease in the Rayleigh number results in a larger influence of the wall-flow. For P r = 5, the number of columnar vortices is smaller (see also the discussions of Figs. 8 and 14), and this is reflected in the melting morphology. For the flow regime considered here, increasing Re also decreases the rate of melting, even as  The parameters are the same as in Fig. 2. Due to the enhanced heat transport, the plots here are at earlier times than shown in Fig. 2.
the Nusselt number increases with increasing Re. This is discussed further in §IV B.

D. Influence of velocity boundary conditions at the lower surface
In §III A - §III C, the fluid 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., [21][22][23][24][25]. Morever, the critical Rayleigh number in Eq. 26 is highest for free-slip top-and bottom-boundaries, and lowest for one free-slip and one no-slip boundary; the case of two no-slip boundaries is intermediate between these cases [16,26]. Despite this, for the parameter ranges considered here, the Nusselt number is higher 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., 21]. Thus, the melting rates are higher when the lower boundary is no-slip, as seen in Fig. 4, which can be compared to the free-slip lower boundary case of Fig. 2. The simulations in §III A- §III D were performed in a container with rigid no-slip lateral boundaries. We examine the influence of the presence of the walls by performing simulations in horizontally periodic domains with the same domain sizes and the same parameters. The results for the case corresponding to that in §III A are shown in Fig. 5. Clearly the wall-mode is absent in periodic geometries and the peripheral streaming flow, present in the confined geometry, is no longer seen here. The columnarvortical flow at P r = 5 and the resultant melting of the solid in a pattern reflecting the presence of these vortices remain unchanged. Note also that the amount of melting is similar in Figs. 1 and 5.

F. Influence of thermal diffusivity in the solid
The thermal diffusivity of the solid governs the amount of heat transported away from the solid-liquid interface and the melt rate (see Eq. 16). In the simulations presented in §III A- §III E,κ s = 1 and the diffusivity of the solid was the same as that of the liquid. Figure 6 shows the effects ofκ s = 1, for two values ofκ s = 0.2 and κ s = 5, with f = 1 (so that the upper boundary is at θ = −1). As expected, the ability of the solid to conduct heat away from the interface affects the rate of melting (see Eq. 16).

A. Flow structures and melting morphology
As we have seen, for P r = 5 and the ranges of Re and Ro considered here, the flow takes the form of columnar vortices, with a peripheral retrograde near-wall current in the confined geometry. These flow structures control the morphology of the melting. For instance, the peripheral streaming current causes increased melting near the walls of the container, and this effect is more prominent when the Rayleigh number is close to the bulk critical Rayleigh number (Eq. 26). This is shown in Fig. 7, where the azimuthally averaged height of the liquid layer is plotted as a function of the distance from the axis of rotation. We note that since the fluid height h(t) is a function of time, so are the Rayleigh number (Eq. 24) and the critical Rayleigh numbers (Eqs. 26  When the Rayleigh numbers are higher, heat transfer occurs primarily through columnar vortices. These vortices are identified by counting the number of isolated regions where The properties of these vortices, which are themselves  functions of the flow parameters, are essential in controlling the melting morphology. Figure 8 shows that the general expectation that the number of columnar vortices increases with Reynolds number and the rotation rate. Of particular relevance to the phase changes, Fig.  9(a) shows that as the number of vortices increases the average area of each vortex decreases. Moreover, this behavior is independent of the flow regimes studied, as evidenced by the parametric collapse onto a single curve. Figure 9(b) shows that, beyond the initial transients, the total vortex area increases linearly with time until eventually the vortices start to merge (see also Fig. 8) and the total vortex area reaches a quasi-steady state. This total vortex area, in the geostrophic convection regime, increases with increasing Re and Ro. These vortices carry heat from the lower boundary to the solid. Therefore, the morphology of the phase bound-   ary (e.g., the average area and number of "holes" or melt void regions melted into the solid) reflects the state of the flow. Indeed, we see in Figs. 10 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.

B. Heat transport and the melting rate
In §II C we noted that the initial and boundary conditions in most of the simulations reported here, except those in §III F, are that the solid is at the melting temperature throughout, viz., θ(t = 0) = 0, and the upper boundary is held at θ = 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 enthalpy conservation, Eq. 21, as where the terms in square brackets are nondimensional. Dividing through by k l ∆T L gives

ReP r St
is the average nondimensional temperature over the simulation volume and 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 Fig. 11. Initially, all the energy supplied to the system from the boundary heats up the liquid. For lower Ro and Re 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 = 25). 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 vigor of convection depends on Ra, Ro and Re, such a balance requires quasi-steady rotating convection. We see in Fig. 11(a) that the quasi-steady state of convection in the fluid described by Eq. 28 breaks down at t = 120 when fluid comes into contact with the upper solid boundary as the solid vanishes. Note further that in our simulations Equation 21 is satisfied exactly, since we work directly with the enthalpy, and the slight mismatch between − ∂θ ∂z z=0 and the sum ReP r Additionally Fig. 11 shows that when the heat transport is sufficiently rapid that the specific heat stored in the convecting fluid is small, there is a nearly steady balance between the heat supplied at the base of the cell and the melt rate. Indeed, as shown in Fig. 12, 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 Eq. 28 using only the first two terms, viz.,

ReP r St
where Nu is the Nusselt number-the total heat flux scaled by the conductive heat flux-across the fluid region.
In §IV A we showed that for most combinations of Re, Ro and P r examined here, the phase boundary is ramified, so that the solid depth varies substantially with horizontal position. In consequence, we see from Fig. 12 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 Section III of 9] when considering the quasisteady balance in Eq. 31. We note, however, that we understand that there are three-dimensional heat fluxes in the interfacial region, which is simpler to treat when the phase boundary has small amplitude variations, such as in the non-rotating case [e.g., 10,27]. Another perspective is that for a vortex-induced highly ramified interface, the interfacial region can be considered as a "mushy layer", as observed in binary systems [3], 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 Fig. 12 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 Eqs. 24, 25 and 26, as where β is in general a function of Ra/Ra bulk c and C is a numerical prefactor that may depend on P r.  [14] found β = 3/4 and [22] found β = 2/7. In the limit of large Ro, that is in the classical nonrotating Rayleigh-Bénard convection regime, one finds, with a different prefactor than in Eq. 32, β = 1/3 up to Ra = 10 15 [28][29][30][31].
When the instantaneous Nusselt number in this quasisteady state obeys Eq. 32, then the solution of Eq. 31 shows that the rate of motion of the solid-liquid interface, h s = H − h, is controlled by the parameter The logarithmic derivatives of Eqs. 32 and 33 show how the net heat flux, Nu, the melt rate and β depend on the flow parameters Re, Ro, P r and St. However, a rigorous quantitative assessment [see e.g., 32] requires ranges of parameters larger than we have available here. We may nonetheless draw some qualitative conclusions by examining the trends of our simulations and thereby provide some guidance for future studies.
The Nusselt number from Eq. 31 is plotted in Fig. 13  (a) for St = 1, wherein, at t ≈ 150, we observe the intuitive suppression of vertical heat flux with rotation rate; here a linear decrease in Nu with decreasing Ro, consis-tent with β = 3/4. Moreover, the intuitive increase of Nu with P r is seen to be approximately linear and thus consistent with β = 3/4 or 3/2. Fig. 13(b) shows Nu for St = 0.05. We see that in a range of the quasi-steady balance, between t = 750 and t = 1000, Nu doubles when Re doubles. Note that although we expect that an increase in inertia will enhance buoyancy fluxes, because the time range over which Nu is constant is also a function of the other flow parameters, we can only exclude the quadratic dependence on Re associated with β = 3 for these data.
Time histories of the volume-averaged height of the solid, h s = L/A − h, the slope of which are given by Eq. 33, as a function of Re, Ro and P r are shown in Fig. 14. Clearly, rotation suppresses average vertical heat transport and hence the melt rate for all β values, but here, despite the increase in Nu shown in Fig. 13(b) we see a very weak dependence on Re. Moreover, in Figs. 14 (a) and (c) for P r = 1, Ra/Ra bulk c is small, and much of the flow is along the lateral boundaries thereby localizing the melt rate and reducing the mean melt rate relative to the P r = 5 case. The dependence on Ro and the weak dependence on P r and Re is consistent with β = 3/4 or 3/2. Finally, in Fig. 15 (c), although the state of the flow is not modified substantially, we see here the strong influence on decreasing the growth rate with an increase in latent heat (here a decrease in St) as is well known in solidification [3].

C. Maximal Phase Boundary Roughness & Maximal Heat Flux
We conclude §IV with the observation that for systems with large latent heat (a small St here), the quasi-steady balance between the heat supplied at the base of the cell and the growth rate that lead to Eq. 31, and was demonstrated in Fig. 13 (b,c), is also associated with the maximally rough phase boundary. For the smallest Stefan numbers studied here, St = 0.05 and 0.1, we see in Figs. 16(a) and (b) that as the maximal quasi-steady Nusselt numbers are reached, the roughness of the solid-liquid interface also reaches a maximal state, as characterized by the standard deviation of the height of the phase boundary as a function of time. For all values of Re the plateau in Nu seen in Figs. 13 (b) and (c) is accompanied by the maximal roughness of the phase boundary. We note that Nu(t) and the maximal interface roughness depend on the parameter η, but the correlation between Nu and σ(h) shown in Figs. 13 and 16 does not. (See also Fig.  22 in Appendix B.) In non-rotating turbulent Rayleigh-Bénard convection, with Dirichlet boundary conditions and periodically rough boundaries, [33,34] showed that, for a given roughness wavelength, there is a ratio of the thermal boundary layer thickness to the roughness amplitude that optimizes Nu. In that situation, the system is in a statistical steady state. Here we find that, in the regime in which Eq. 31 is valid, as Nu reaches a maximum, the phase boundary geometry evolves towards the roughest state, after which it coarsens while maintaining a quasisteady state flux balance. Whereas the solid void and vortex area (Fig. 10) as well as the magnitude of Nu (Figs. 13 (b) and (c)) depend on the flow parameters, for large latent heat systems the maximum in the roughness is insensitive to those parameters. Thus, although the roughness becomes maximal when the maximal Nu is reached, the maximal Nu state persists as the interfacial region coarsens. Here, as opposed to the case where the roughness is static, the phase boundary is free to evolve to a geometry determined by the underlying conservation laws. However, whether this correlation is associated with a variational principle, wherein the interfacial geometry evolves to maximize the heat transport, or the geometry is simply that which pins heat transporting vortices, cannot presently be distinguished.

V. CONCLUSIONS
We have studied the melting of a pure solid from its liquid phase when the former overlies the latter, which is heated from below and the entire system is rotated about an axis parallel to gravity. The width of the system is twice its depth and we have examined ranges of the  is given by Eq. 27 [19]. Third, in the periodic geometry, there is no flow for Ra < Ra bulk c . We found that the number of melt voids in the solid is proportional to the num-ber of heat transporting vortices present, which in turn increases with increasing Reynolds number and decreasing Rossby number. However, the overall melting rate is a nontrivial function of the flow parameters. For example, although the number of vortices and the melt regions they produce increase with Reynolds number, the average area of each vortex decreases with Reynolds number and thus so too does the average melt rate. Generally the average melt rate increases with increasing Rossby and Prandtl numbers and with decreasing Reynolds number. Moreover, 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. 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, characterized by the Nusselt number, 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., 2, 3, 35, 36]. Moreover, it has been shown experimentally, mathematically and numerically that boundary roughness enhances heat transport in otherwise classical non-rotating turbulent Rayleigh-Bénard convection with Dirichlet boundary conditions [34,37,38]. Here, in a pure system, we find that convective and rotationally controlled vortices alone can create ramified phase boundaries. Moreover, for systems with a large latent heat of fusion, the maximal heat transport and the maximal phase boundary roughness occur at the same time. We suggest that in astrophysical and geophysical problems wherein rotational effects are important, it may be prudent to modify analytical treatments of the associated Stefan problems when planarity of the phase boundary is assumed. acknowledged.
Computations were performed on Tetralith. The Swedish Research Council under grant no. 638-2013-9243, is gratefully acknowledged for support.

APPENDIX A : VALIDATION OF THE ENTHALPY METHOD
The enthalpy method used here has been validated by comparing it to the one-dimensional analytical solution for a purely conducting case [e.g., 3]. Consider a semiinfinite solid layer in the region z > 0 at the melting temperature. The boundary at z = 0 is held at θ = 1. The solid melts, forming a liquid layer of height h(t) given by where ξ is the solution of the transcendental equation deriving from the Stefan condition, In Fig. 17 (a) the analytical solution of the Stefan problem is compared with a numerical solution of Eq. 21 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 θ = 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 θ = 0. The boundaries are held at θ(z = 0) = 1 and θ(z = H) = 0. The numerical solution in one dimension is compared with the solution from the 3D solver, and the amount of unmelted solid plotted as a function of time in Fig. 17(b). In both these cases, the 1D solution is obtained using fourth-order Runge-Kutta integration; the 3D solver uses a second-order Adams-Bashforth scheme (as described in §II D). 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 < χ < 1, and is distributed over a finite number of gridpoints. This is shown in Fig. 18(a) where the enthalpy φ is plotted on a vertical line through the peak an inverted hole in the solid region. The mask function χ varies from 0 to 1 over a distance of about δz = 0.015, which is 8 gridpoints in the 512 2 × 256 simulations. We note that the thinness of the interface is not affected by changes in the grid-resolution or in the penalization parameter, as seen from Fig. 18(b), with η = 5 × 10 −4 , and a coarser resolution.

APPENDIX B: PENALIZATION PARAMETER
The volume penalization method has a tunable parameter η. All of our results are reported with the penalization parameter η = 10 −3 ( §II D). 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 η creates a velocity boundary layer of size (νη) 1/2 in the solid. [39] showed that the optimal value of η is such that the grid spacing is comparable to the boundary layer thickness, namely dx ∼ (νη) 1/2 .
In detail the melting process is influenced by the boundary layer and hence depends on η. However, as seen in Fig. 19(a), upon reduction of η by a factor of 10, the melt rate changes by only a few percent. Therefore, the latent heat flux and the quasi-steady balance described by Eq. 31 underlying the results shown in Figs. 11(b), 13 and 16 are insensitive to the choice of η.
Finally, whereas some qualitative features of the phase boundary arise when η is changed, snapshots of the interface shown in Fig. 20 demonstrate the persistence of the central behavior; convective vortices etch voids into the solid, and the number of voids are proportional to the number of vortices (Fig. 21). Thus, as noted in §IV C, Nu(t) and the maximal interface roughness depend on η, but the correlation between Nu and σ(h) shown in Figs. 13, 16 and 22 do not.    Fig. 19(c,d). (Cf., Figs. 13, 16.)