1. Introduction
Fluid motions that accompany solid–liquid phase transitions are ubiquitous in both the natural and engineering environments (Epstein & Cheung Reference Epstein and Cheung1983; Glicksman, Coriell & McFadden Reference Glicksman, Coriell and McFadden1986; Huppert Reference Huppert1986; Worster Reference Worster2000; Meakin & Jamtveit Reference Meakin and Jamtveit2009; Hewitt Reference Hewitt2020). Such fluid motions either are a result of the generation of unstable density gradients during solidification (Davis, Müller & Dietsche Reference Davis, Müller and Dietsche1984; Dietsche & Müller Reference Dietsche and Müller1985; Wettlaufer, Worster & Huppert Reference Wettlaufer, Worster and Huppert1997; Worster Reference Worster1997; Davies Wykes et al. Reference Davies Wykes, Huang, Hajjar and Ristroph2018), or are externally imposed and thereby influence morphological and/or hydrodynamic instabilities in the system (Delves Reference Delves1968, Reference Delves1971; Gilpin, Hirata & Cheng Reference Gilpin, Hirata and Cheng1980; Coriell et al. Reference Coriell, McFadden, Boisvert and Sekerka1984; Forth & Wheeler Reference Forth and Wheeler1989; Feltham & Worster Reference Feltham and Worster1999; Neufeld & Wettlaufer Reference Neufeld and Wettlaufer2008a,Reference Neufeld and Wettlauferb; Dallaston, Hewitt & Wells Reference Dallaston, Hewitt and Wells2015; Ramudu et al. Reference Ramudu, Hirsh, Olson and Gnanadesikan2016; Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019). The interactions between fluid flows and phase-changing surfaces underlie our understanding of, among other phenomena, the evolution of Earth's cryosphere (McPhee Reference McPhee2008; Hewitt Reference Hewitt2020), defects in materials castings (Kurz & Fischer Reference Kurz and Fischer1984), and production of single crystals for silicon chips (Worster Reference Worster2000). In this study, we focus on understanding the combined effects of buoyancy, rotation and mean shear on the evolution of a phase boundary of a pure material.
The directional solidification of a pure melt is qualitatively different from that of a multi-component melt. However, there are qualitative similarities in the responses of the two systems to externally imposed shear flows (Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019). During the solidification of a multi-component melt, solute particles (atoms, molecules or colloids) are rejected into the bulk liquid, which can result in the constitutional supercooling of the liquid adjacent to the phase boundary (Worster, Peppin & Wettlaufer Reference Worster, Peppin and Wettlaufer2021). This promotes ‘morphological instability’, which leads to dendritic growth: the formation of a convoluted crystal matrix with a concentrated solution of solute particles trapped in the interstices (Worster Reference Worster2000). This region is called a mushy layer. Depending on the equation of state of the liquid phase and the growth direction relative to the gravitational field, convection within (the ‘mushy layer mode’) or adjacent to (the ‘boundary layer mode’) the mushy layer may result (Worster Reference Worster1992). Externally imposed flows can change the nature of these two modes.
Some of the earliest systematic investigations into the effects of a mean shear flow on the directional solidification of binary alloys are due to Delves (Reference Delves1968, Reference Delves1971) and Coriell et al. (Reference Coriell, McFadden, Boisvert and Sekerka1984). Delves (Reference Delves1968, Reference Delves1971) considered the effects of a parabolic flow on morphological instability using linear stability analysis in two dimensions. He found that the parabolic flow suppresses the instability, the degree of which depends on material properties. Similarly, Coriell et al. (Reference Coriell, McFadden, Boisvert and Sekerka1984) studied the effects of a Couette flow on the morphological and thermo-solutal instabilities in three dimensions, and found that the latter are suppressed more than the former, which is due to the fact that the wavelength of morphological instability is several orders of magnitude smaller than the wavelength of thermo-solutal instability. Furthermore, the use of a Couette flow as the base-state velocity profile in this case may be valid only locally, since the momentum equations are governed by an asymptotic suction boundary layer profile (Drazin & Reid Reference Drazin and Reid2004; Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019).
Forth & Wheeler (Reference Forth and Wheeler1989) studied the effects of an asymptotic suction boundary layer flow on morphological instability during the directional solidification of a binary alloy using linear stability analysis in three dimensions. They found that under certain conditions, the shear flow inhibits the morphological instability and may lead to the development of travelling waves along the phase boundary. However, the phase boundary was found to have negligible effects on the stability of the shear flow itself.
In the directional solidification of mushy layers, the fate of any incipient perturbation introduced at the mush–melt interface depends naturally on the interaction between the mushy and boundary layer modes of convective motion. In order to understand the evolution of such perturbations, Feltham & Worster (Reference Feltham and Worster1999, see also the corrigendum in Neufeld et al. Reference Neufeld, Wettlaufer, Feltham and Worster2006) considered the effects of flows of inviscid and viscous melts on the mush–melt interface using linear stability analysis in two dimensions. Neglecting buoyancy, they found that the forced flows over a corrugated interface introduce pressure perturbations at the phase boundary, which in turn drive fluid motions in the mushy layer. These fluid motions, in certain cases, were found to amplify the interfacial perturbations. However, the heat flux from the melt was found to have negligible effects on this amplification, and was found to generate only travelling waves at the interface.
Neufeld & Wettlaufer (Reference Neufeld and Wettlaufer2008a,Reference Neufeld and Wettlauferb) used a combination of linear stability analysis in three dimensions and experiments to examine the effects of a shear flow on the convective motions in the mushy layer and bulk melt. Their main findings are: (i) below a critical strength of the imposed shear flow, the convective motions are moderately suppressed; (ii) above this critical strength, the stability of these convective motions decreases monotonically with increasing shear flow strength; and (iii) for a sufficiently strong shear flow, quasi-two-dimensional striations of zero volume fraction of solid are formed perpendicular to the shear flow direction at the mush–melt interface. These striations form due to local dissolution and growth of the mushy layer, which result from the interactions between the shear and convective motions within and outside the mushy layer.
In addition to shear flows, several studies have explored the effects of rotation on the stability of the mush–melt interface and the convective motions inside and outside the mushy layer. For example, Sample & Hellawell (Reference Sample and Hellawell1982, Reference Sample and Hellawell1984) explored experimentally the effects of rotation and rotational-and-precessional motions during directional solidification of 
${\rm NH}_4{\rm Cl}\text {--}{\rm H}_2{\rm O}$ and Pb–Sn systems, respectively. Their key findings are as follows. (i) Channels in the mushy layer through which solute is expelled into the bulk melt originate close to the mush–melt interface due to perturbations in the solutal boundary layer. (ii) When the mould is rotated by approximately an angle to the vertical of 
$20^{\circ }$–
$30^{\circ }$ at low speeds (
${<}5$ rpm), the induced interfacial shear flow arrests the growth of the perturbations, thereby inhibiting the formation of channels. (iii) Rotation about the vertical axis has very little effect on the evolution of the perturbations. In a related study, Lu & Chen (Reference Lu and Chen1997) used linear stability analysis to show that the rotation about a vertical axis has an appreciable stabilizing effect only for very large rotation rates (
${\sim }10^5$ rpm).
Motivated by the experiments of Sample & Hellawell (Reference Sample and Hellawell1984), Chung & Chen used linear stability analysis to study the detailed effects of inclined rotation and precessional motions on the stability of the mushy layer (Chung & Chen Reference Chung and Chen2000), the mush–melt interface and the bulk melt (Chung & Chen Reference Chung and Chen2003). They found the following. (i) Rotation of the system about an inclined axis induces flow parallel to the mush–melt interface in both the mushy layer and the bulk melt. However, the maximum velocity in the bulk is much larger than the maximum velocity in the mushy layer (Chung & Chen Reference Chung and Chen2000). (ii) The induced velocity in the melt increases with increasing angle of inclination and decreases with rotation rate, but the induced velocity in the mushy layer is sensitive only to the inclination angle and increases with it (Chung & Chen Reference Chung and Chen2000). (iii) The stabilization of the mushy layer is due largely to the reduction in buoyancy perpendicular to the mush–melt interface (Chung & Chen Reference Chung and Chen2000). (iv) When there is only precessional motion, the most unstable mode of instability in the mushy layer is in the direction perpendicular to the mush–melt interface. Introducing rotation allows the induced flow to explore all directions equally, thereby stabilizing all modes equally (Chung & Chen Reference Chung and Chen2000). (v) Five competing mechanisms – the (stabilizing) reduction of buoyancy and rotation normal to the interface, the (destabilizing) component of gravity parallel to the interface, and the (stabilizing or destabilizing) flow induced by inclination and precession – determine the stability of the bulk melt and the interface (Chung & Chen Reference Chung and Chen2003).
There have been relatively few studies on the effects of buoyancy, mean shear and rotation on directional solidification of pure melts. Using linear and weakly nonlinear stability analyses and experiments, Davis et al. (Reference Davis, Müller and Dietsche1984) explored the different patterns on the phase boundary that result due to thermal convection. These patterns and the transition between them were further explored experimentally by Dietsche & Müller (Reference Dietsche and Müller1985). Recent studies have focused on the effects of high-Rayleigh-number convection on the phase boundary (Esfahani et al. Reference Esfahani, Hirata, Berti and Calzavarini2018; Favier, Purseed & Duchemin Reference Favier, Purseed and Duchemin2019; Purseed et al. Reference Purseed, Favier, Duchemin and Hester2020). A detailed summary of these studies can be found in Toppaladoddi (Reference Toppaladoddi2021).
Hirata, Gilpin & Cheng (Reference Hirata, Gilpin and Cheng1979a), Hirata et al. (Reference Hirata, Gilpin, Cheng and Gates1979b) and Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) studied experimentally the effects of laminar and turbulent boundary layer flows on the evolution of phase boundaries between ice and water. Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) observed that when the conductive heat flux through the ice is less than approximately twice the water-to-ice flux, a perturbation introduced initially at the interface grew and advected downstream, leading to the formation of ‘rippled’ surfaces. The heat transfer rate increased by approximately 30–60 % over these corrugated surfaces when compared with a planar interface. These observations were attributed solely to mean shear by Gilpin et al. (Reference Gilpin, Hirata and Cheng1980), but a recent re-analysis of their data showed that there were substantial buoyancy effects due to anomalous behaviour of water at 4 
$^\circ$C in their experiments (Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019). Furthermore, a linear stability analysis of the Rayleigh–Bénard–Couette flow over a phase boundary showed that mean shear has a stabilizing effect and buoyancy has a destabilizing effect on the phase boundary (Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019). More recently, the interplay between mean shear, buoyancy and phase boundaries has been explored further in two dimensions (Toppaladoddi Reference Toppaladoddi2021) and three dimensions (Couston et al. Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2020) using direct numerical simulations.
Ravichandran & Wettlaufer (Reference Ravichandran and Wettlaufer2021) explored the combined effects of rotation and thermal convection on the evolution of a phase boundary in three dimensions. They found that the rotation-dominated convection of a pure liquid, which takes the form of columnar vortices (Rossby Reference Rossby1969), melts voids into its overlying solid phase. They showed that the feedback of the melting on the convective flow can arrest the horizontal motion of flow structures like the peripheral streaming flow in wall-bounded rotating convection, and hypothesized that the columnar vortices may also be pinned to the locations of the voids. Furthermore, their study revealed that the degree of interfacial roughness co-varies with the heat flux.
Here, we use direct numerical simulations to study the melting of the solid phase of a pure substance driven by thermal convection of its melt, subject to a constant applied pressure gradient in the horizontal direction and rotation about the vertical axis. The combined effects of shear, rotation and buoyancy on the evolution of phase boundaries are important in several geophysical settings (Wells & Wettlaufer Reference Wells and Wettlaufer2008; Ramudu et al. Reference Ramudu, Hirsh, Olson and Gnanadesikan2016; Hewitt Reference Hewitt2020), and provide important motivations for this study.
The remainder of the paper is organized as follows. In § 2, we formulate the problem and present the equations of motion along with the initial and boundary conditions. We discuss briefly the numerical method used to solve the equations of motion in § 3. We then discuss results from simulations of non-rotating flows in § 4.1, and of rotation-influenced flows in § 4.2. Finally, we conclude with a summary of our observations and suggestions for future work in § 5.
2. Problem formulation and governing equations
In this study, we consider a cuboid of dimensions 
$L \times L \times 2H$. The aspect ratio of the cell is defined as 
$A = L/2H$ and is fixed at either 
$4$ or 
$8$ in the simulations reported here. Figure 1 shows a vertical cross-section of the domain. At the initial time, the planar phase boundary is located at 
$z = h_0$ and the fluid is contained in the region 
$0 \le z \le h_0$. By definition, the temperature at the phase boundary is 
$T = T_m$, and that at the top plate, 
$z = 2H$, is maintained at 
$T = T_m$. The bottom plate located at 
$z = 0$ is maintained at a temperature above the melting point, 
$T = T_m + \Delta T$. We note that because of the isothermal temperature field in the solid phase, it completely melts before a stationary state is reached. Thus we focus solely on the system dynamics until the solid phase vanishes. Periodic boundary conditions are imposed in the horizontal directions.
Figure 1. A vertical cross-section of the geometry considered. The size of the domain is 
$L \times L\times 2H$, with aspect ratio 
$A=L/2H$. The bottom surface is maintained at a higher temperature than the interface. No-slip and no-penetration conditions for the fluid velocity are imposed at the bottom boundary and the evolving interface, and periodic boundary conditions are imposed in the horizontal.
The equations of motion in the different regions of the domain are as discussed below.
2.1. Liquid
The conservations of momentum, mass and heat are given by
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\frac{1}{\rho_0} \, \boldsymbol{\nabla} p + g\alpha(T-T_m)\, \hat{\boldsymbol{z}} - 2\varOmega\,\hat{\boldsymbol{z}} \times \boldsymbol{u} + \nu \, {\nabla}^2 \boldsymbol{u} + F_p \, \hat{\boldsymbol{x}}, \end{gather}
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \text{{and}}\end{gather}
\begin{gather} \frac{\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \kappa \, {\nabla}^2T, \end{gather}
respectively. Here, 
$\boldsymbol {u}(\boldsymbol {x},t) = (u, v, w)$ is the three-dimensional velocity field; 
$\rho _0$ is the reference density; 
$p(\boldsymbol {x},t)$ is the pressure field; 
$g$ is acceleration due to gravity; 
$\alpha$ is the thermal expansion coefficient; 
$T(\boldsymbol {x},t)$ is the temperature field; the constant pressure gradient is applied as a horizontal body force per unit mass, 
$F_p$; 
$\boldsymbol {\hat {x}}$ and 
$\boldsymbol {\hat {z}}$ are the unit vectors along the 
$x$- and 
$z$-axes, respectively; 
$\varOmega$ is the constant speed of rotation; 
$\nu$ is the kinematic viscosity; and 
$\kappa$ is the thermal diffusivity.
2.2. Solid
The entire solid phase is isothermal because of the temperature boundary condition at the upper boundary. Hence the temperature field in the solid phase is given by
\begin{equation} T(\boldsymbol{x},t) = T_m. \end{equation}2.3. Evolution of the phase boundary
The location of the phase boundary is determined using the Stefan condition, which is a statement of energy balance (Worster Reference Worster2000):
\begin{equation} \rho_0 \, L_s \, v_n ={-} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{q_l},\quad \text{at} \ z = h(x,y,t). \end{equation}
Here, 
$L_s$ is the latent heat of fusion; 
$v_n$ is the normal component of growth rate of the solid phase; 
$\boldsymbol {n}$ is the unit normal pointing into the liquid; and 
$\boldsymbol {q_l}$ is the heat flux towards the phase boundary from the liquid.
2.4. Boundary conditions
We impose Dirichlet conditions on the temperature at the bottom and top boundaries of the domain:
\begin{equation} T(x, y, z=0,t) = T_m + \Delta T \quad \text{and} \quad T(x,y,z = H,t) = T_m, \end{equation}
where 
$\Delta T > 0$. The temperature at the phase boundary is the bulk melting temperature:
\begin{equation} T(x,y,z=h,t) = T_m. \end{equation}We impose no-slip and no-penetration conditions on the flow velocity at the lower boundary and the phase boundary, given by
\begin{gather} u(x,y,z = 0, t) = v(x,y,z=0,t) = w(x,y,z = 0, t) = 0, \end{gather}
\begin{gather} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{t}_1 = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{t}_2 = 0,\quad\text{at} \ z = h(x,y,t), \end{gather}
respectively, where 
$\boldsymbol {t}_1$ and 
$\boldsymbol {t}_2$ are the 
$x$ and 
$y$ projections of the unit tangents at the phase boundary. For simplicity, we assume that the reference densities of the solid and liquid phases are equal to 
$\rho _0$. Hence there is no additional normal velocity created at the phase boundary due to density difference between the phases. We also impose periodic boundary conditions for the temperature and velocity fields at 
$x = 0$ and 
$x = L$, and 
$y = 0$ and 
$y=L$.
2.5. Initial conditions
At 
$t=0$, the temperature of the liquid layer is maintained at 
$T(\boldsymbol {x},t=0) = T_m$, and the velocity field is 
$\boldsymbol {u}(\boldsymbol {x},t=0) = \boldsymbol {0}$. The shear flow in the horizontal direction is driven by the body force starting at 
$t=0$.
2.6. Non-dimensional equations of motion
To non-dimensionalize the equations of motion and the associated boundary conditions, we choose 
$H$ as the length scale, 
$U_0 = \sqrt {g\alpha \, \Delta T \, H}$ as the velocity scale, 
$t_0 = H/U_0$ as the time scale, 
$p_0 = \rho _0U_0^2$ as the pressure scale, and 
$\Delta T$ as the temperature scale. Using these in (2.1)–(2.3) and (2.5), and maintaining pre-scaled notation, we have
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla} p + \sqrt{\frac{Pr}{Ra}} \, {\nabla}^2 \boldsymbol{u} - \frac{2}{Ro} \, \hat{\boldsymbol{z}} \times \boldsymbol{u} + \theta \, \hat{\boldsymbol{z}} + {\frac{1}{Ri_b}} \, \hat{\boldsymbol{x}}, \end{gather}
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather} \frac{\partial \theta}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta = \sqrt{\frac{1}{Pr \times Ra}} \, {\nabla}^2 \theta \ , \text{ and} \end{gather}
\begin{gather} v_n = \sqrt{\frac{1}{Pr \times Ra \times St^2}} \, \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta,\quad \text{at} \ z = h(x,y,t), \end{gather}where the dimensionless temperature,
\begin{equation} \theta = \frac{T-T_m}{\Delta T}, \end{equation}
is defined such that 
$\theta =0$ at the phase boundary, and 
$\theta =1$ at the heated lower boundary. The five dimensionless parameters that govern the dynamics of the system are the Rayleigh (
$Ra$), bulk Richardson (
$Ri_b$), Prandtl (
$Pr$), Stefan (
$St$) and Rossby (
$Ro$) numbers, defined as
\begin{gather} Ra = \frac{g\alpha \, \Delta T \, H^3}{\nu\kappa}, \quad Ri_b = \frac{g\alpha \, \Delta T}{F_p}, \quad Pr = \frac{\nu}{\kappa}, \end{gather}
\begin{gather} St = \frac{L_s}{C_p \, \Delta T} \quad \text{and} \quad Ro = \sqrt{\frac{g\alpha \, \Delta T}{\varOmega^2H}}, \end{gather}
where 
$C_p$ is the specific heat of the solid phase. These represent the ratios of buoyancy to viscous forces (
$Ra$), of buoyancy to the imposed horizontal body force per unit mass (
$Ri_b$), of momentum to heat diffusivities (
$Pr$), of the latent to specific heats (
$St$), and of rotational to convective time scales (
$Ro$). Note that 
$St$ is defined based on the temperature difference across the liquid layer and not the solid layer. A bulk Reynolds number can be defined by combining 
$Ra$ and 
$Ri_b$, and is given by
\begin{equation} Re = \sqrt{\frac{Ra}{Ri_b \times Pr}}. \end{equation}
In cases where both rotation and shear are present, we define 
$\varSigma$ as the ratio of centrifugal to applied body forces,
\begin{equation} \varSigma = \frac{Ri_b}{Ro^2} = \frac{\varOmega^2 H}{F_p}, \end{equation}
such that mean shear dominates when 
$\varSigma < 1$.
The dimensionless temperature boundary conditions are
\begin{equation} \theta(x, y, z=0,t) = 1 \quad \text{and} \quad \theta(x,y,z = H,t) = 0 \end{equation}at the bottom and top surfaces, respectively, and
\begin{equation} \theta(x,y,z=h,t) = 0 \end{equation}at the phase boundary. The no-slip and no-penetration boundary conditions for the velocity field remain unaltered after the non-dimensionalization.
An important means of quantifying the response of the system is in terms of the Nusselt number, which is defined as the ratio of total heat flux to the heat flux due only to conduction. We define the Nusselt number as
\begin{equation} Nu = {\sqrt{Ra\,Pr\,St^2} \, h \, \frac{{\rm d}h}{{\rm d}t}}, \end{equation}
where 
$h$ is the horizontally averaged height of the fluid layer. It is the total heat flux scaled by the conductive heat flux over the liquid layer, whose depth is evolving. This follows Ravichandran & Wettlaufer (Reference Ravichandran and Wettlaufer2021), although their Stefan number is the inverse of that defined here.
3. Numerical method
Equations (2.11)–(2.13), along with the boundary conditions, are solved using the finite-volume solver Megha-5, with a volume penalization method for the melting solid. The solver uses second-order central differences in space and a second-order Adams–Bashforth time-stepping scheme, and has been validated extensively for free shear flows (Singhal et al. Reference Singhal, Ravichandran, Diwan and Brown2021, Reference Singhal, Ravichandran, Govindarajan and Diwan2022) and Rayleigh–Bénard convection (Ravichandran, Meiburg & Govindarajan Reference Ravichandran, Meiburg and Govindarajan2020; Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2020, Reference Ravichandran and Wettlaufer2021). Details of the volume-penalization algorithm – including validation against the analytical solution for the purely conductive Stefan problem, tests of sensitivity to the penalization parameter, and the convergence of the solution under grid-refinement – may be found in Ravichandran & Wettlaufer (Reference Ravichandran and Wettlaufer2021). To test the results for grid independence in this study, simulations for the lowest 
$Ri_b$ cases were repeated with double the number of grid points in the vertical, and the difference in the results was found to be negligible. Furthermore, we have verified that the body force prescribed in (2.1) produces the known analytical plane-Poiseuille velocity profile with an exponential relaxation time scale of 
$({\rm \pi} ^2 \nu )^{-1}$ for the peak velocity.
4. Results
Here, we present our results on the effects of mean shear, buoyancy and rotation on the evolution of the phase boundary. A list of the parameter combinations studied is presented in Tables 1. In order to highlight the effects of rotation, we first consider only the effects of mean shear and buoyancy, and later introduce rotation. Choosing 
$Pr > 1$ ensures that 
$Ro < 1$ without unduly increasing the resolution requirement for the Ekman boundary layers at the upper and lower surfaces. For this reason, we choose 
$Pr = 5$ for the cases where rotation is present (see e.g. Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2020). In the absence of rotation, there is no qualitative difference in the melting dynamics for 
$Pr = 1$ and 
$5$. The choice of 
$Ra$ in this study reflects a balance between choosing a sufficiently large 
$Ra$ to observe the convective dynamics whilst minimizing the computational effort required to resolve the thermal boundary layers.
Table 1. List of cases reported here. In all cases, the Rayleigh number is 
$Ra=1.25\times 10^5$, the Stefan number is 
$St=1$, and the penalization parameter is 
$\eta =1.4\times 10^{-3}$. We also give the bulk Reynolds number 
$Re$ for the non-rotating cases (2.17), and the parameter 
$\varSigma$ for the rotating cases (2.18). The solutions obtained are independent of the grid resolution, time step 
$\Delta t$, and volume penalization parameter 
$\eta$.
4.1. Zero rotation (
$Ro=\infty$)
When rotational effects are absent, the evolution of the phase boundary is affected by the combined action of the mean shear flow and convection. The relative strength of mean shear to buoyancy is quantified using 
$Ri_b$: buoyancy dominates when 
$Ri_b \gg 1$, and mean shear dominates when 
$Ri_b \ll 1$. We explore the range 
$Ri_b \in [1, \infty )$ in this section, with the other dimensionless parameters fixed at 
$Ra=1.25\times 10^5$ and 
$Pr = St = 1$.
In figures 2(a–d), we show the effects of increasing the strength of the mean shear flow on the phase boundary at 
$t=85$.
Figure 2. The solid–liquid interface viewed from the solid side for 
$Ro=\infty$, 
$Ra=1.25\times 10^5$, 
$Pr=1$ and 
$St=1$ at 
$t=85$ for (a) 
$Ri_b = \infty$, (b) 
$Ri_b=100$, (c) 
$Ri_{b}=10$, (d) 
$Ri_{b}=1$. The arrow in (a) denotes the direction of the shear flow along the 
$x$-direction. For 
$Ri_b=\infty$, the melting creates dome-like voids in the solid, as seen in (a). Nascent streamwise alignment of the melting pattern appears for 
$Ri_{b} = 100$ in (b). For smaller 
$Ri_b$, the formation of quasi-two-dimensional grooves separated by sharp furrows is seen in (c,d). The mean liquid layer depths are comparable in (a–c), while the liquid depth in (d) is approximately 
$20\,\%$ larger (see figure 7).
When mean shear is absent (
$Ri_b = \infty$), convective motions tend to create dome-like features at the boundary, which ‘lock in’ the convection cells. This is seen in figure 2(a). Such features have been investigated in detail in both two (Favier et al. Reference Favier, Purseed and Duchemin2019; Purseed et al. Reference Purseed, Favier, Duchemin and Hester2020) and three (Esfahani et al. Reference Esfahani, Hirata, Berti and Calzavarini2018) dimensions. As the strength of the mean shear is increased (decreasing 
$Ri_b$), quasi-two-dimensional corrugations with a wavelength transverse to the direction of the shear flow begin to emerge at the phase boundary, as seen in figures 2(b–d). The shear flow, which is parallel to the 
$x$-axis, tends to homogenize the phase boundary along the flow direction. However, it does not affect the corrugations perpendicular to the flow, parallel to the 
$y$-axis (Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019). Such features have also been observed in the recent study of Couston et al. (Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2020). The effect of the mean shear flow here is consistent with that in two dimensions, where the corrugations of a phase boundary decrease in amplitude as a result of shear (Toppaladoddi Reference Toppaladoddi2021).
The quasi-two-dimensional features at the phase boundary are due to the emergence of the longitudinal rolls, which are the preferred form of convection when 
$Ra > Ra_c$ and 
$Re$ is not too large (Clever & Busse Reference Clever and Busse1991). These longitudinal rolls can be discerned in the contours of the vertical components of the velocity and temperature. While these effects begin to appear for 
$Ri_b=100$ (not shown), they are clearly evident for 
$Ri_b=10$, as shown in figure 3. Vertical sections of the flow in figures 3(c,d) show that the interface is homogeneous in the direction of the shear flow. These phase boundary patterns are qualitatively different from those observed experimentally by Gilpin et al. (Reference Gilpin, Hirata and Cheng1980). This is due to the fact that in their experiments, 
$Re$ based on the boundary layer thickness was 
$O(10^4)$, whereas in our simulations bulk 
$Re$ is 
$O(10^2\text {--}10^3)$, and the wavelength of transverse patterns is 
$O(10^3)$ times the viscous sublayer thickness (Claudin, Durán & Andreotti Reference Claudin, Durán and Andreotti2017) and hence longer than our simulation domain. Longitudinal features are also observed in flows over eroding (Allen Reference Allen1971) and dissolving (Cohen et al. Reference Cohen, Berhanu, Derr and Pont2016, Reference Cohen, Berhanu, Derr and Pont2020) boundaries, where the associated grooves eventually merge to produce scallops (see also Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019). Another interesting feature to note is that for 
$Ri_b = 10$, the corrugations at the phase boundary become more regular in the 
$y$–
$z$ plane, which is shown in figure 3( f). This shows that the mean shear flow inhibits vertical motions, which is similar to its effects in two dimensions (Toppaladoddi Reference Toppaladoddi2021).
Figure 3. Contours of vertical velocity 
$w$ and temperature 
$\theta$ at the horizontal section 
$z=H/2$ (a,b) and the vertical sections 
$y=0$ (c,d) and 
$x=0$ (e, f) at 
$t=84$. The solid–liquid interface is shown in the vertical sections as a solid yellow line. The flow parameters are 
$Ri_{b}=10$, 
$Ro=\infty$, 
$Ra=1.25\times 10^5$, 
$Pr=1$ and 
$St=1$. Here, quasi-two-dimensional corrugations with a wavelength transverse to the direction of the shear flow are clearly associated with longitudinal roll structures in the flow.
The longitudinal rolls persist down to 
$Ri_b \approx 4$, below which they start to merge. The effects of the convective rolls on the melting morphology are shown in the Hövmöller diagrams of the liquid layer depth as functions of 
$y$ and 
$t$ in figure 4. At 
$t=150$, for 
$Ri_b = 10$, the six corrugations have wavelength approximately 2.66, which is slightly larger than twice the initial fluid depth, whereas for 
$Ri_b=4$, the merging and coarsening of corrugations is evident (cf. figures 4a,b). Prior to this time, the dynamics of the corrugations depends on 
$Ri_b$, highlighting the importance of the specific history in controlling the interfacial structure. The coupling of the flow structures and the melting morphology is thus seen more starkly here than in convection without shear in non-rotating (e.g. Esfahani et al. Reference Esfahani, Hirata, Berti and Calzavarini2018) and rotating (Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2021) cases, where there are domes, like those in figure 2(a), rather than longitudinal structures, whose average size increases continuously as the solid melts.
Figure 4. Hövmöller plots of the fluid height at 
$x=0$, plotted as a function of 
$(t,y)$ for (a) 
$Ri_b=10$, showing the emergence of six domes of nearly equal widths (of approximately 
$2.66$), and (b) 
$Ri_b=4$, showing the emergence of three domes each of two distinct sizes. In both instances, the domes maintain their horizontal positions.
When 
$Ri_b$ is decreased further, the longitudinal rolls lose coherence and the bulk flow becomes turbulent and fully three-dimensional. This is shown in figures 5 and 6(a,b) for 
$Ri_b = 1$. Nonetheless, remnants of the longitudinal rolls can still be seen along the phase boundary parallel to the 
$y$-axis (figures 5e, f), which is homogeneous along the 
$x$-axis (figure 5c). Clearly, as seen in figure 6(c), the classical log layer of a turbulent shear flow is responsible for a significant temperature gradient in the bulk. In contrast, for 
$Ri_b=10$, the bulk is nearly isothermal as expected in convection-dominated flows.
Figure 5. Contours of vertical velocity 
$w$ and temperature 
$\theta$ as in figure 3 and for the same flow parameters except 
$Ri_{b}=1$. The increased shear, compared to 
$Ri_b=10$, causes the flow to become fully three-dimensional. This greatly increases the heat transfer rate and the Nusselt number (see also figure 7). As a result, the grooves are less prominent (cf. figure 2).
Figure 6. Profiles of the horizontally averaged (a) horizontal velocity 
$\bar {u}$ and (b) temperature 
$\bar {\theta }$, at 
$t\approx 85$ for 
$Ri_b=1$ and 
$Ri_b=10$ , with 
$Ro=\infty$, 
$Ra=1.25\times 10^5$, 
$Pr=1$, 
$St=1$. The velocity and thermal boundary layer thicknesses (
$30$ wall units for the velocity) are shown using dashed lines; the average height of the fluid layer is shown using dotted lines; the boundary layers are thinner for smaller 
$Ri_b$. (c) For sufficiently small 
$Ri_b$, a turbulent boundary layer develops at the lower (smooth) boundary, with the well-known logarithmic velocity profile for 
$z^+ \gtrsim 30$, where 
$z^+ = z/\delta _\tau$ is the vertical coordinate scaled in terms of the wall-friction thickness 
$\delta _\tau = \nu / u_\tau$, and 
$u_\tau = [\nu (\partial \bar {u} / \partial z )_{z=0} ]^{1/2}$ is the friction velocity. As a result, for large 
$Ri_b$, convection dominates, and the fluid bulk is nearly isothermal; for small 
$Ri_b$, on the other hand, the fluid bulk has a constant temperature gradient. The friction Reynolds numbers 
$Re_\tau = u_\tau H / \nu$, based on the initial liquid height, are 
$85.4$ and 
$306.5$ for 
$Ri_b=10$ and 
$Ri_b=1$, respectively.
The melting rates and interface roughness are functions of the flow properties and, as we have seen from figure 4, can also determine the behaviour of flow structures. In figure 7(a) we show the time evolution of 
$h_s$ and 
$Nu$ for different 
$Ri_b$. The melting rate of the solid, ascertained from the slopes of the 
$h_s(t)$ curves, is a non-monotonic function of 
$Ri_b$. When 
$Ri_b$ is reduced from 
$\infty$ to 
$4$, the rate of melting decreases. However, with a further decrease in 
$Ri_b$ from 
$4$ to 
$1$, there is an increase in the melting rate. This non-monotonicity results from the relative dominance of shear flow and convection in the dynamics of the system, and can also be seen in figures 7 and 8. For 
$Ri_b = 100$, the mean shear flow is relatively weak, but still acts to inhibit vertical motions. This results in decreased heat transport relative to the 
$Ri_b = \infty$ case. This trend continues for 
$Ri_b = 10$, where the vertical motions are further inhibited. However, a further decrease in 
$Ri_b$ to 
$1$ results in the mean shear flow becoming dominant, and 
$Nu$ and the melting rate increase due to the action of forced convection. These effects of decreasing 
$Ri_b$ can be seen by comparing figures 3( f) and 5( f), where the flow field transitions from having distinct plume structures to a more chaotic three-dimensional flow. The effects are qualitatively consistent with those observed in the study of pure Rayleigh–Bénard–Poiseuille flow in three dimensions by Scagliarini, Gylfason & Toschi (Reference Scagliarini, Gylfason and Toschi2014).
Figure 7. (a) Average thickness of the solid layer 
$h_s$ (red dashed curves, left axis) and the melting Nusselt number (from (2.21); blue curves, right axis). (b) Roughness of the solid–liquid interface 
$\sigma _h$ for 
$Ro=\infty$, 
$Ra=1.25\times 10^5$, 
$Pr=1$, 
$St=1$, and varying 
$Ri_b$. The rate of melting, and thus the Nusselt number, varies non-monotonically with 
$Ri_b$, as seen in (a). Similarly, the interface roughness also varies non-monotonically with 
$Ri_b$. (See also figure 8.)
The roughness of the solid–liquid interface, defined here as the standard deviation 
$\sigma _h$ of the fluid height 
$h(x,y)$, is plotted in figure 7(b). Similar to the melting rate, we see that the roughness also varies non-monotonically with increasing 
$Ri_b$, and is largest in the absence of shear, 
$Ri_b=\infty$. For 
$Ri_b = 100$, the incipient homogenization of the flow structures, and thus the melting morphology, causes a decrease in 
$\sigma _h$. The emergence of the streamwise corrugations commensurate with convective rolls offsets the decrease of the roughness. As a result, an increase in the maximum of 
$\sigma _h$ is seen at 
$Ri_b=10$. For 
$Ri_b\leq 10$, the dominance of shear and the resulting turbulent flow lead to less prominent morphological features and hence a reduction in 
$\sigma _h$. Furthermore, we note that by comparing figures 7(a,b), it can be seen that the roughness and the Nusselt number reach their maximal values as the solid–liquid interface reaches the upper boundary, as seen in the absence of shear (Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2021).
4.2. The effects of rotation
Having studied the effects of mean shear and buoyancy on the flow dynamics, we now consider the effects of rotation. We use the same 
$Ra=1.25\times 10^5$ as in § 4.1, and choose 
$Pr=5$ to give 
$Ro = 0.6325$, which is fixed for all the cases discussed in this subsection. The range of 
$Ri_b$ explored here is 
$Ri_b \in [0.1, 100]$, which is slightly larger than in § 4.1.
Figures 9(a–d) show the phase boundary geometries for the different 
$Ri_b$ when the solid layer has nearly completed melted. For 
$Ri_b = 100$, the effects of shear are small and the phase boundary has the regular dome-like features seen in purely rotating case, where columnar vortices span the height of the liquid layer and transport heat between the lower and upper boundaries (Rossby Reference Rossby1969; Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2021), as seen in figure 9(a). Upon decreasing 
$Ri_b$, the homogeneous distribution of domes is replaced by a weak transverse structure, as seen for 
$Ri_b = 10$ in figure 9(b). With a further decrease in 
$Ri_b$, the phase boundary develops patterns that are similar to the ones seen in absence of rotation (compare figures 2c,d and 9c,d). The alignment of the longitudinal corrugations/grooves differs between rotating and non-rotating cases because in the former case the flow in the bulk is determined by the balance between the Coriolis effect and the applied pressure gradient. Therefore, the mean shear flow in the bulk is in the negative 
$y$-direction rather than in the positive 
$x$-direction. Hence, despite the external pressure gradient applied in the positive 
$x$-direction, the longitudinal grooves in the presence of rotation are parallel to the 
$y$-axis.
Figure 8. Maximum values of (a) 
$Nu$ and (b) 
$\sigma _h$ in figure 7 as functions of 
$Ri_b$.
Figure 9. The solid–liquid interface seen from the solid side with flow parameters 
$Ro = 0.6325$, 
$Ra= 1.25\times 10^5$, 
$Pr=5$, 
$St=1$, and (a) 
$Ri_{b}=100$ at 
$t=155$, (b) 
$Ri_{b}=10$ at 
$t=155$, (c) 
$Ri_{b}=1$ at 
$t=240$, (d) 
$Ri_{b}=0.5$ at 
$t=240$. (a) Characteristic melting by rotating convection (see e.g. Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2021) is seen for large 
$Ri_b$. Compare the effect of large 
$Ri_b$ for 
$Ro=\infty$ in figure 2(b), and note that the figures are plotted at different times owing to the differences in the time scales of melting. (b) For 
$Ri_{b}=10$, the lateral drift of the vortices is reflected in the connectivity between neighbouring solid voids to form longer grooves that do not span the width of the domain. (c) For 
$Ri_{b}=1$, the grooves become quasi-two-dimensional and span the width of the domain in the lateral direction. (d) For 
$Ri_{b}=0.5$, the grooves are less prominent than in (c). The times chosen are such that the liquid layer depths are comparable. The values of 
$\varSigma$ for these cases are (a) 250, (b) 25, (c) 2.5, and (d) 1.25, respectively. Hence, in all these cases, rotation dominates over mean shear.
The merging of the columnar vortices for 
$Ri_b = 1$ can also be seen in the vertical velocity and temperature fields, as shown in figure 10. In particular, figure 10 panels (a) and (b) show some of the merger events and the incipient longitudinal streaks that develop as a result of the mean shear flow. These features can also be seen in the vertical cross-sections of the flow (figure 10c–f).
Figure 10. Contours of vertical velocity 
$w$ and temperature 
$\theta$ at the horizontal section 
$z=H/2$ (a,b), and the vertical sections 
$y=0$ (c,d) and 
$x=0$ (e, f). The solid–liquid interface is shown in the vertical sections as a solid yellow line. The flow parameters are 
$Ri_b=1$, 
$Ro = 0.6325$, 
$Ra=1.25\times 10^5$, 
$Pr=5$ and 
$St=1$. The value of 
$\varSigma$ is 
$0.25$ for this case, so that mean shear dominates rotation. The drift velocity of the columnar vortices along the negative 
$y$-direction is sufficiently strong for the melting morphology to become quasi-two-dimensional (compare the vertical sections (e, f) to the melting morphology shown in figure 9c).
The loss of coherence of the columnar vortices is complete for 
$Ri_b \le 0.2$, and the longitudinal rolls now appear clearly in the flow field. This is seen in figure 11, which shows the contours of vertical velocity and temperature for 
$Ri_b = 0.1$. In addition to the complete merger of the vortices, one can also see that the longitudinal rolls meander about their mean location. This indicates the onset of a wave-like instability of the rolls that is qualitatively similar to that seen in thermal convection with mean shear flows, but without rotation (Clever & Busse Reference Clever and Busse1991, Reference Clever and Busse1992; Pabiou, Mergui & Bénard Reference Pabiou, Mergui and Bénard2005). These waves are clearly visible in the patterns at the phase boundary, as seen in figure 12(d).
Figure 11. Contours of vertical velocity 
$w$ and temperature 
$\theta$ as in figure 10, and for the same parameters except 
$Ri_b=0.1$ (
$\varSigma =0.25$) and initial liquid layer height 
$h_0=0.5$. The columnar vortices are no longer distinguishable and have merged to produce two-dimensional rolls in the bulk (see also figure 13). The turning of the flow in the Ekman layers creates the grooves aligned at an angle intermediate between the 
$x$- and 
$y$-directions.
Figure 12. Morphology of the solid–liquid interface (seen from the solid side) at (a) 
$t=354$ and (b) 
$t=495$, for 
$Ri_{b}=0.2$, and (c) 
$t=495$ and (d) 
$t=552$, for 
$Ri_{b}=0.1$. Here, 
$\varSigma < 1$ for both cases, and mean shear also dominates rotation. As melting proceeds, the grooves are aligned in a direction intermediate between the 
$x$- and 
$y$-directions, with a smaller angle between the grooves and the 
$x$-axis for smaller 
$Ri_{b}$. At late times, the interface develops a sinusoidal mode parallel to the grooves. The initial liquid height is 
$h_{0}=0.5$ in these simulations.
To understand the alignment of the longitudinal corrugations with respect to the applied pressure gradient, we study the mean flow in the bulk. In figure 13, we plot the horizontal components of the mean velocity averaged over the horizontal planes and the angle that the mean horizontal velocity makes with the applied pressure gradient for 
$Ri_b=0.1$. These plots show that despite the strength of the applied pressure gradient, the flow in the bulk remains in geostrophic balance (thus the angle made by the flow relative to the 
$x$-direction is 
$\phi =-{\rm \pi} /2$). However, the no-slip upper and lower boundaries force the mean flow direction to rotate and form Ekman spirals. Thus this counter-clockwise rotation of the mean flow direction (
$\phi$ increases from 
$-{\rm \pi} /2$ to 
$0$) in the Ekman layer adjacent to the upper boundary has profound effects on the melting morphology. This is seen in figure 12, which shows the formation of grooves at the phase boundary for 
$Ri_b = 0.2$ (figures 12a,b) and 
$Ri_b = 0.1$ (figures 12c,d). In contrast to the large 
$Ri_b$ cases, these grooves are aligned neither parallel nor perpendicular to the direction of the applied pressure gradient, but at an intermediate angle.
Figure 13. (a) Horizontally averaged mean values of the velocities 
$u$ (left axis, blue) and 
$v$ (right axis, red) along the 
$x$- and 
$y$-directions, respectively. (b) Angle 
$\phi =\tan ^{-1}(\bar {v}/\bar {u})$ relative to the applied pressure gradient for 
$Ri_b=0.1$ at 
$t=552$ corresponding to the snapshot in figure 12(d). The fluid velocity is nearly perpendicular to the applied pressure gradient in the bulk, and turns counter-clockwise in the Ekman layer at the solid–liquid interface. For 
$Ri_b/Ro^2<1$, this counter-clockwise turning of the flow in the upper Ekman layer leads to the observed formation of grooves at an oblique angle. Note that the grid point at 
$z=0$ is not shown.
The Ekman layer at the upper boundary may be defined as the region where 
$\phi \neq -{\rm \pi} /2$ (see figure 13b). For smaller 
$Ri_b$, the average horizontal velocity in this region is larger and makes a slightly smaller angle with the 
$x$-axis, decreasing from 
$85^\circ$ for 
$Ri_b=0.2$, to 
$80^\circ$ for 
$Ri_b=0.1$ (at the times shown in figure 12). Similarly, we find that the angle between the grooves and the 
$x$-axis decreases slightly, from 
$75^\circ$ to 
$70^\circ$, as 
$Ri_b$ decreases.
In the absence of rotation (figure 2), and with finite background rotation for 
$Ri_b \geq {O}(1)$ (figure 9), the grooves formed are stationary in that once formed they do not move, but become more prominent with time (see figure 4). In contrast, the obliquely aligned grooves for small 
$Ri_b$ are advected perpendicular to their lengths, as seen from the Hövmöller plots in figure 14 for 
$Ri_b=0.1$. Since the propagation speeds are proportional to the wavelengths in the 
$x$- and 
$y$-directions, figures 14(a,b) also show that the grooves translate horizontally perpendicular to their lengths.
Figure 14. Hövmöller plots of the fluid height at (a) 
$x=0$ and (b) 
$y=0$, showing the drift of the grooves for the same parameters as in figure 11. The grooves move perpendicular to their lengths and thus maintain the same angle to the 
$x$-axis for the duration of their existence.
The impact of the flow and phase boundary dynamics on the time evolution of 
$h_s$, 
$Nu$ and 
$\sigma _h$ for the cases corresponding to figure 9 are shown in figure 15. The melt rate and the Nusselt number decrease as 
$Ri_b$ decreases, because as the strength of the shear flow increases, vertical convective motions are suppressed, as shown in figure 15(a). For large 
$Ri_b$, columnar convection dominates the flow, leading to a more ramified dome structure of the phase boundary, which gives way to an aligned periodic phase boundary as 
$Ri_b$ decreases and shear dominates. Thus 
$\sigma _h$ increases with 
$Ri_b$, as shown in figure 15(b). Note that for 
$Ri_b = 10$ and 
$100$, 
$\sigma _h$ reaches a maximum at the same time as the corresponding 
$Nu$, as observed in the case of 
$Ri_b = \infty$ (Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2021).
Figure 15. (a) Thickness of the solid layer 
$h_s$ (red dashed curves, left axis) and the melting Nusselt number (from (2.21); blue curves, right axis). (b) Roughness of the solid–liquid interface 
$\sigma _h$ for 
$Ri_b=100,10,1,0.5$ (with 
$\varSigma = 250, 25, 2.5, 1.25$, respectively) and other parameters as in figure 10. The rate of melting decreases with decreasing 
$Ri_{b}$, unlike in figure 7. We also note that the maxima of 
$Nu$ and 
$\sigma _h$ are reached simultaneously for 
$Ri_b=10$.
Similar behaviour in the evolution of 
$h_s$, 
$Nu$ and 
$\sigma _h$ is observed for the smallest three 
$Ri_b$ cases studied, as shown in figure 16. However, for large times, there are substantial jumps in 
$Nu$ and 
$\sigma _h$ for 
$Ri_b = 0.1$ and 
$0.2$. This is due to the onset of the wave-like interfacial instability of the longitudinal rolls shown in figures 11 and 12. This change in the phase boundary geometry leads to the generation of more intense downwelling plumes that result in enhanced heat transfer (Toppaladoddi, Succi & Wettlaufer Reference Toppaladoddi, Succi and Wettlaufer2015). Another effect of these cold plumes is that the bulk fluid is at a temperature lower than the mean of the top and bottom surfaces, as shown in figure 17 for 
$Ri_b = 0.1$.
Figure 16. (a) The thickness of the solid layer 
$h_s$ and the Nusselt number Nu and (b) the roughness of the phase boundary 
$\sigma_h$ as a function of time, for 
$Ro = 0.6325$, 
$Ra=1.25\times 10^5$, 
$Pr=5$, 
$St=1$, with 
$\varSigma =1.25, 0.5, 0.25$ for 
$Ri_b=0.5,0.2,0.1$, respectively. The overall rate of melting is initially smaller for smaller 
$Ri_{b}$, but as the flow develops, the Nusselt number increases sharply with the development of the along-groove sinusoidal mode. Also note the abrupt change in the slope of the thickness curve for 
$Ri_b=0.1$ at 
$t \approx 552$, which coincides with the increase in 
$Nu$.
Figure 17. Profiles of the horizontally averaged (a) horizontal velocity 
$-\bar {v}$ and (b) temperature 
$\bar {\theta }$, for 
$Ro=0.6325$, 
$Ra=1.25\times 10^5$, 
$Pr=5$, 
$St=1$, 
$Ri_b=0.1$ at the times shown. The temperature in the fluid bulk has a pronounced temperature gradient due to the combined influence of shear and rotation in the initial stages of evolution; the emergence of the along-groove sinusoidal mode coincides with the disappearance of the bulk temperature gradient. The dashed lines show the thermal boundary layer at the lower surface, while the dotted lines show the average height of the fluid layer.
In geostrophic convection, the interaction between the thermal and momentum boundary layers plays an important role in the transport of heat (Rossby Reference Rossby1969; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012; King, Stellmach & Aurnou Reference King, Stellmach and Aurnou2012). In figure 17, we plot the vertical profiles of the horizontally averaged temperature 
$\bar {\theta }$ and the horizontally averaged velocity 
$-\bar {v}$ for 
$Ri_b=0.1$. At early times, the thermal boundary layer is much thicker than the Ekman layer and the heat transport is controlled by rotation. This also results in a significant bulk temperature gradient. As the depth of the liquid layer increases with time, the effective Rayleigh number increases and the thermal boundary layer becomes thinner, and is comparable to the Ekman layer thickness at 
$t\approx 500$. The crossing of the thermal and velocity boundary layers leads to convection-dominated dynamics, which is evident from the small temperature gradient in the bulk at 
$t\approx 550$ shown in figure 17(b).
The different phase boundary geometries suggest the existence of four distinct regimes. (1) For 
$Ri_b \gg 1$, the mean shear flow is weak, and rotation and buoyancy dominate, leading to columnar vortices and dome-like features at the phase boundary. (2) For 
$Ri_b = O(1)$, the mean shear flow is stronger, and the columnar vortices are advected with the flow. This is reflected in the phase boundary geometry, which shows incipient grooves perpendicular to the direction of the applied pressure gradient (figures 9b,c and 10a,b). (3) For 
$Ri_b \lesssim 0.5$, strong shear leads to the loss of coherence of the columnar vortices, resulting in their partial merger, while the interfacial grooves become less prominent. This is seen in figure 9(d). (4) For 
$Ri_b \lesssim 0.2$, the columnar vortices merge completely, and convective motions become quasi-two-dimensional in the bulk. The flow in these cells is turned counter-clockwise in the Ekman layers, leading to the formation of grooves aligned at an angle to the direction of the bulk shear flow and thereby become sinusoidal.
5. Conclusions
We have systematically explored the effects of buoyancy, rotation and shear on the evolution of a phase boundary using direct numerical simulations in three dimensions for 
$Ro = \infty\ {\rm and}\ 1$, 
$Ra=1.25\times 10^5$, 
$St = 1$, 
$Pr = 1\ {\rm and}\ 5$, and 
$Ri_b \in [0.1, \infty )$. The main conclusions from our study are as follows.
(1) In the absence of rotation (
$Ro = \infty$), we observe either dome-like features or grooves on the phase boundary depending on whether buoyancy or mean shear dominates the flow. In particular, we find the following three features.
(i) As the value of
$Ri_b$ decreases from $\infty$, the strength of the shear flow increases, and both the flow structure and phase boundary geometry transition from being three-dimensional to being quasi-two-dimensional. This is seen in figures 2 and 3(a,b).(ii) For small
$Ri_b$, the grooves formed on the phase boundary are aligned parallel to the direction of the shear flow. These grooves are due to longitudinal rolls, which are the preferred form of convection in this case (Clever & Busse Reference Clever and Busse1991). This is in agreement with the numerical results of Couston et al. (Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2020), but in contrast to the experimental results of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980). In our simulations, the bulk Reynolds number is $Re = O(10^2\text {--}10^3)$, whereas in the experiments of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980), $Re$ based on the boundary layer thickness is $O(10^4)$, so the flow regimes are not directly comparable.(iii) The dimensionless heat flux,
$Nu$, is a non-monotonic function of $Ri_b$. For $Ri_b \in [10, 100]$, the mean shear is weaker than buoyancy, but nonetheless inhibits vertical motions, leading to smaller melt rates and vertical heat transport. On further decreasing $Ri_b$ to $1$, the flow becomes three-dimensional and enters the forced turbulent convection regime. This results in larger heat transport and hence an increased melt rate.
(2) When rotation is introduced (
$Ro=0.6325$), the mean flow in the interior is in geostrophic balance and thus in a direction perpendicular to the applied pressure gradient. The dynamics are governed by the balance between mean shear and rotation, and are represented by the parameter $\varSigma$ (see (2.18)), with the imposed shear becoming dominant for $\varSigma <1$. The following features are observed for the rotation-influenced flows.
(i) For
$Ri_b = 100$ ($\varSigma \gg 1$), the shear flow is weak and the phase boundary consists of dome-like features.(ii) On reducing
$Ri_b$ to $10$ ($\varSigma = 25$), these dome-like features begin to merge. A further reduction in $Ri_b$ leads to the formation of grooves at the phase boundary.(iii) For
$Ri_b \in [0.2, 0.5]$ ($\varSigma \in [ 0.5, 1.25 ]$), the merger of the dome-like features is complete. The rolls develop a wave-like instability for $Ri_b \!=\! 0.2$ and $0.1$ ($\varSigma \!=\! 0.5$ and $0.25$, respectively), which leads to sinusoidal evolution of the phase boundary. This instability is similar qualitatively to that observed in thermal convection in the presence of a mean shear flow (Clever & Busse Reference Clever and Busse1991, Reference Clever and Busse1992; Pabiou et al. Reference Pabiou, Mergui and Bénard2005).
Our study provides foundational understanding of the effects of buoyancy, mean shear and rotation on phase-changing boundaries. Although we have treated values of 
$Ra$ and 
$Re$ that are much larger than previously examined, the parameter range explored here is representative of many but certainly not all geophysical settings (see e.g. Meakin & Jamtveit Reference Meakin and Jamtveit2009; Neufeld, Goldstein & Worster Reference Neufeld, Goldstein and Worster2010; Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019; Weady et al. Reference Weady, Tong, Zidovska and Ristroph2022, and references therein). Our results suggest several avenues for further research. With or without rotation, our findings suggest that a secondary instability disrupts the basic flow state in which the applied pressure gradient is balanced by other forces in the flow. The mechanism of the instability, and the boundary between the quasi-two-dimensional state and the three-dimensional state, are compelling foci for further study. In particular, in rotating convection, the competition between the geostrophic bulk and the viscous Ekman layers is crucial in determining the state of the flow and the concomitant heat transport, which is a question of significant interest in the literature (King et al. Reference King, Stellmach and Aurnou2012; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012). Here, however, this competition is reflected in the morphology of the solid–liquid interface, with either the bulk flow or the Ekman boundary layers controlling the alignment of the observed grooves. Therefore, the manner in which the free boundary reflects the state of the flow, and the associated flow-geometry feedback, provide an ideal testing ground for general questions in rotating convection and a framework for the wide range of settings, from geophysical to industrial, in which such flows accompany phase changes.
Funding
Computational resources from the Swedish National Infrastructure for Computing (SNIC) under grants SNIC/2019-3-386 and SNIC/2020-5-471 are gratefully acknowledged. Computations were performed on Tetralith. S.R. and J.S.W. acknowledge support from Swedish Research Council under grant no. 638-2013-9243, and S.T. acknowledges a Research Fellowship from All Souls College, Oxford. Nordita is partially supported by Nordforsk.
Declaration of interests
The authors report no conflict of interest.
Open access
