* Kazemi and Poole [2018] Seyed Mehran Kazemi and David Poole.  Simple embedding for link prediction in knowledge graphs.  In Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada, pages 4289–4300, 2018.
  * Lacroix et al. [2018] Timothée Lacroix, Nicolas Usunier, and Guillaume Obozinski.  Canonical tensor decomposition for knowledge base completion.  In Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, volume 80 of Proceedings of Machine Learning Research, pages 2869–2878. PMLR, 2018.
  * Marino et al. [2017] Kenneth Marino, Ruslan Salakhutdinov, and Abhinav Gupta.  The more you know: Using knowledge graphs for image classification.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, pages 20–28. IEEE Computer Society, 2017.
  * Meilicke et al. [2019] Christian Meilicke, Melisachew Wudage Chekol, Daniel Ruffinelli, and Heiner Stuckenschmidt.  Anytime bottom-up rule learning for knowledge graph completion.  In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI 2019, Macao, China, August 10-16, 2019, pages 3137–3143. ijcai.org, 2019.
  * Mondal et al. [2021] Sutapa Mondal, Sumit Bhatia, and Raghava Mutharaju.  EmEL++: Embeddings for EL++ description logic.  In Proceedings of the AAAI 2021 Spring Symposium on Combining Machine Learning and Knowledge Engineering (AAAI-MAKE 2021), Stanford University, Palo Alto, California, USA, March 22-24, 2021, volume 2846 of CEUR Workshop Proceedings. CEUR-WS.org, 2021.
  * Omran et al. [2018] Pouya Ghiasnezhad Omran, Kewen Wang, and Zhe Wang.  Scalable rule learning via learning representation.  In Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, IJCAI 2018, July 13-19, 2018, Stockholm, Sweden, pages 2149–2155. ijcai.org, 2018.
  * Pavlovic and Sallinger [2023] Aleksandar Pavlovic and Emanuel Sallinger.  ExpressivE: A spatio-functional embedding for knowledge graph completion.  In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net, 2023.
  * Reinanda et al. [2020] Ridho Reinanda, Edgar Meij, and Maarten de Rijke.  Knowledge graphs: An information retrieval perspective.  Found. Trends Inf. Retr., 14(4):289–444, 2020.
  * Schneider et al. [2022] Phillip Schneider, Tim Schopf, Juraj Vladika, Michael Galkin, Elena Simperl, and Florian Matthes.  A decade of knowledge graphs in natural language processing: A survey.  In Proceedings of the 2nd Conference of the Asia-Pacific Chapter of the Association for Computational Linguistics and the 12th International Joint Conference on Natural Language Processing, AACL/IJCNLP 2022 - Volume 1: Long Papers, Online Only, November 20-23, 2022, pages 601–614. Association for Computational Linguistics, 2022.
  * Sun et al. [2019] Zhiqing Sun, Zhi-Hong Deng, Jian-Yun Nie, and Jian Tang.  Rotate: Knowledge graph embedding by relational rotation in complex space.  In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net, 2019.
  * Toutanova and Chen [2015] Kristina Toutanova and Danqi Chen.  Observed versus latent features for knowledge base and text inference.  In Proceedings of the 3rd Workshop on Continuous Vector Space Models and their Compositionality, CVSC 2015, Beijing, China, July 26-31, 2015, pages 57–66. Association for Computational Linguistics, 2015.
  * Trouillon et al. [2017] Théo Trouillon, Christopher R. Dance, Éric Gaussier, Johannes Welbl, Sebastian Riedel, and Guillaume Bouchard.  Knowledge graph completion via complex tensor factorization.  J. Mach. Learn. Res., 18:130:1–130:38, 2017.
  * Trouillon et al. [2019] Théo Trouillon, Éric Gaussier, Christopher R. Dance, and Guillaume Bouchard.  On inductive abilities of latent factor models for relational learning.  J. Artif. Intell. Res., 64:21–53, 2019.
  * Wang et al. [2018] Yanjie Wang, Rainer Gemulla, and Hui Li.  On multi-relational link prediction with bilinear models.  In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), the 30th innovative Applications of Artificial Intelligence (IAAI-18), and the 8th AAAI Symposium on Educational Advances in Artificial Intelligence (EAAI-18), New Orleans, Louisiana, USA, February 2-7, 2018, pages 4227–4234. AAAI Press, 2018.
  * Wang et al. [2023] Dingmin Wang, Yeyuan Chen, and Bernardo Cuenca Grau.  Efficient embeddings of logical variables for query answering over incomplete knowledge graphs.  In Thirty-Seventh AAAI Conference on Artificial Intelligence, AAAI 2023, Thirty-Fifth Conference on Innovative Applications of Artificial Intelligence, IAAI 2023, Thirteenth Symposium on Educational Advances in Artificial Intelligence, EAAI 2023, Washington, DC, USA, February 7-14, 2023, pages 4652–4659. AAAI Press, 2023.
  * Xiong et al. [2022] Bo Xiong, Nico Potyka, Trung-Kien Tran, Mojtaba Nayyeri, and Steffen Staab.  Faithful embeddings for EL++ knowledge bases.  In The Semantic Web - ISWC 2022 - 21st International Semantic Web Conference, Virtual Event, October 23-27, 2022, Proceedings, volume 13489 of Lecture Notes in Computer Science, pages 22–38. Springer, 2022.
  * Zhang et al. [2019] Shuai Zhang, Yi Tay, Lina Yao, and Qi Liu.  Quaternion knowledge graph embeddings.  In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pages 2731–2741, 2019.
  * Zhu et al. [2021] Zhaocheng Zhu, Zuobai Zhang, Louis-Pascal A. C. Xhonneux, and Jian Tang.  Neural bellman-ford networks: A general graph neural network framework for link prediction.  In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, pages 29476–29490, 2021.


# The combined effects of buoyancy, rotation, and shear on phase boundary
evolution

S. Ravichandran  Nordita, KTH Royal Institute of Technology and Stockholm
University, Stockholm 10691, Sweden  S. Toppaladoddi  University of Oxford,
Oxford OX1 4AL  J. S. Wettlaufer  Nordita, KTH Royal Institute of Technology
and Stockholm University, Stockholm 10691, Sweden  Yale University, New Haven,
Connecticut 06520-8109, USA

###### Abstract

We use well-resolved numerical simulations to study the combined effects of
buoyancy, pressure-driven shear and rotation on the melt rate and morphology
of a layer of pure solid overlying its liquid phase in three dimensions at a
Rayleigh number $Ra=1.25\times 10^{5}$. During thermal convection we find that
the rate of melting of the solid phase varies non-monotonically with the
strength of the imposed shear flow. In the absence of rotation, depending on
whether buoyancy or shear dominates the flow, we observe either domes or
ridges aligned in the direction of the shear flow respectively. Furthermore,
we show that the geometry of the phase boundary has important effects on the
magnitude and evolution of the heat flux in the liquid layer. In the presence
of rotation, the strength of which is characterized by the Rossby number,
$Ro$, we observe that for $Ro={O}(1)$ the mean flow in the interior is
perpendicular to the direction of the constant horizontal applied pressure
gradient. As the magnitude of this pressure gradient increases, the geometry
of solid-liquid interface evolves from the voids characteristic of melting by
rotating convection, to grooves oriented perpendicular or obliquely to the
direction of the pressure gradient.

##  I Introduction

Fluid motions that accompany solid-liquid phase transitions are ubiquitous in
both the natural and engineering environments (Epstein and Cheung, 1983;
Glicksman  _et al._ , 1986; Huppert, 1986; Worster, 2000; Meakin and Jamtveit,
2009; Hewitt, 2020). Such fluid motions are either a result of the generation
of unstable density gradients during solidification (Davis  _et al._ , 1984;
Dietsche and Müller, 1985; Wettlaufer  _et al._ , 1997; Worster, 1997; Davies
Wykes  _et al._ , 2018), or are externally imposed and thereby influence
morphological and/or hydrodynamic instabilities in the system (Delves, 1968,
1971; Gilpin  _et al._ , 1980; Coriell  _et al._ , 1984; Forth and Wheeler,
1989; Feltham and Worster, 1999; Neufeld and Wettlaufer, 2008a, b; Dallaston
_et al._ , 2015; Ramudu  _et al._ , 2016; Bushuk  _et al._ , 2019). The
interactions between fluid flows and phase-changing surfaces underlie our
understanding of, among other phenomena, the evolution of Earth’s cryosphere
(McPhee, 2008; Hewitt, 2020), defects in materials castings (Kurz and Fischer,
1984), and production of single crystals for silicon chips (Worster, 2000). 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 response of the two systems to externally imposed shear flows
(Toppaladoddi and Wettlaufer, 2019). 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  _et al._ , 2021). 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, 2000). 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, 1992). 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 (1968, 1971) and Coriell  _et al._ (1984). Delves (1968, 1971)
considered the effects of a parabolic flow on morphological instability using
linear stability analysis in two dimensions (2D). He found that the parabolic
flow suppresses the instability, the degree of which depended on material
properties. Similarly, Coriell  _et al._ (1984) studied the effects of a
Couette flow on the morphological and thermo-solutal instabilities in three
dimensions (3D), 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 only be locally valid, since the
momentum equations are governed by an asymptotic-suction boundary-layer
profile (Drazin and Reid, 2004; Toppaladoddi and Wettlaufer, 2019).

Forth and Wheeler (1989) 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 3D. 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 naturally depends on the
interaction between the mushy– and boundary–layer modes of convective motion.
In order to understand the evolution of such perturbations, Feltham and
Worster (1999, see also the corrigendum in ()) considered the effects of flows
of inviscid and viscous melts on the mush-melt interface using linear
stability analysis in 2D. 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 only generate
travelling waves at the interface.

Neufeld and Wettlaufer (2008a, b) used a combination of linear stability
analysis in 3D 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 of the mushy layer. For example, Sample and
Hellawell (1982, 1984) experimentally explored the effects of rotation and
rotational-and-precessional motions during directional solidification of
NH4Cl-H2O and Pb-Sn systems, respectively. Their key findings are: (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 rotates about an angle to the
vertical ($None$ \- $None$) at low speeds ($<5$ rpm), the induced interfacial
shear flow arrests the growth of the perturbations, thereby inhibiting the
formation of channels; and (iii) rotation about the vertical axis has very
little effect on the evolution of the perturbations. In a related study, Lu
and Chen (1997) 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 and Hellawell (1984), 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 and Chen,
2000), the mush-melt interface and the bulk melt (Chung and Chen, 2003). They
found that (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 and Chen, 2000); (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 only sensitive
to the inclination angle and increases with it (Chung and Chen, 2000); (iii)
the stabilization of the mushy layer is largely due to the reduction in
buoyancy perpendicular to the mush-melt interface (Chung and Chen, 2000); (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 equally stabilizing all modes (Chung and Chen,
2000); (v) five competing mechanisms – the (stabilising) reduction of buoyancy
and rotation normal to the interface, the (destabilising) component of gravity
parallel to the interface, and the (stabilising or destabilising) flow induced
by inclination and precession – determine the stability of the bulk melt and
the interface (Chung and Chen, 2003).

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._ (1984)
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 and Müller (1985). Recent studies
have focussed on the effects of high Rayleigh-number convection on the phase
boundary (Esfahani  _et al._ , 2018; Favier  _et al._ , 2019; Purseed  _et
al._ , 2020). A detailed summary of these studies can be found in Toppaladoddi
(2021).

Hirata  _et al._ (1979a, b) and Gilpin  _et al._ (1980) experimentally studied
the effects of laminar and turbulent boundary-layer flows on the evolution of
phase boundaries between ice and water. Gilpin  _et al._ (1980) observed that
when the conductive heat flux through the ice is less than approximately twice
the water-to-ice flux, an initially introduced perturbation at the interface
grew and advected downstream, leading to the formation of “rippled” surfaces.
The heat transfer rate increased by about 30-60% over these corrugated
surfaces when compared with a planar interface. These observations were
attributed solely to mean shear by Gilpin  _et al._ (1980), but a recent
reanalysis of their data showed that there were substantial buoyancy effects
due to anomalous behaviour of water at 4 ∘C in their experiments (Toppaladoddi
and Wettlaufer, 2019). 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 and Wettlaufer, 2019). More recently, the interplay
between mean shear, buoyancy, and phase boundaries has been further explored
in 2D (Toppaladoddi, 2021) and 3D (Couston  _et al._ , 2020) using direct
numerical simulations.

Ravichandran and Wettlaufer (2021) explored the combined effects of rotation
and thermal convection on the evolution of a phase boundary in 3D. They found
that the rotation-dominated convection of a pure liquid, which takes the form
of columnar vortices (Rossby, 1969), 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 and Wettlaufer, 2008; Ramudu  _et al._ , 2016; Hewitt, 2020),
and provide important motivations for this study.

The remainder of the paper is organised as follows. In §II, we formulate the
problem and present the equations of motion along with the initial and
boundary conditions. We briefly discuss the numerical method used to solve the
equations of motion in §III. We then discuss results from simulations of non-
rotating flows in §IV.1 and of rotation-influenced flows in §IV.2. Finally, we
conclude with a summary of our observations and suggestions for future work in
§V.

##  II Problem formulation and governing equations

In this study, we consider a cuboid of dimensions of $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.

Figure 1:  A vertical cross-section of the geometry considered. The size of
the domain is $L\times L\times 2H$, with the 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.

At the initial instant, the planar phase boundary is located at $z=h_{0}$ and
the fluid is contained in the region $0\leq z\leq 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.

The equations of motion in the different regions of the domain are as follows.

###  II.1 Liquid

The conservation of momentum, mass and heat are given by

$\frac{\partial\boldsymbol{u}}{\partial
t}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}=-\frac{1}{\rho_{0}}\,\nabla
p+g\,\alpha\,\left(T-T_{m}\right)\,\boldsymbol{\hat{z}}-2\,\Omega\boldsymbol{\hat{z}}\times\boldsymbol{u}+\nu\,\nabla^{2}\boldsymbol{u}+F_{p}\,\boldsymbol{\hat{x}},$
(1) $\nabla\cdot\boldsymbol{u}=0,\text{{and}}$ (2) $\frac{\partial T}{\partial
t}+\boldsymbol{u}\cdot\nabla T=\kappa\,\nabla^{2}T,$ (3)

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; $\Omega$ is the constant speed of rotation; $\nu$ is the
kinematic viscosity; and $\kappa$ is the thermal diffusivity.

###  II.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

$T(\boldsymbol{x},t)=T_{m}.$ (4)

###  II.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, 2000):

$\rho_{0}\,L_{s}\,v_{n}=-\boldsymbol{n}\cdot\boldsymbol{q_{l}}\hskip
5.69046pt\text{at}\hskip 2.84544ptz=h(x,y,t).$ (5)

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.

###  II.4 Boundary conditions

We impose Dirichlet conditions on temperature at the bottom and top boundaries
of the domain:

$T(x,y,z=0,t)=T_{m}+\Delta T\quad\mbox{and}\quad T(x,y,z=H,t)=T_{m},$ (6)

where $\Delta T>0$. The temperature at the phase boundary is the bulk melting
temperature:

$T(x,y,z=h,t)=T_{m}.$ (7)

We impose no-slip and no-penetration conditions on the flow velocity at the
lower boundary and the phase boundary, given by

$u(x,y,z=0,t)=v(x,y,z=0,t)=w(x,y,z=0,t)=0;$ (8)
$\boldsymbol{u}\cdot\boldsymbol{n}=\boldsymbol{u}\cdot\boldsymbol{t_{1}}=\boldsymbol{u}\cdot\boldsymbol{t_{2}}=0\quad\mbox{at}\quad
z=h(x,y,t),$ (9)

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
induced 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$.

###  II.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$.

###  II.6 Non-dimensional equations of motion

To non-dimensionalise 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_{0}\,U_{0}^{2}$ and the pressure scale; and $\Delta T$ as the
temperature scale. Using these in equations 1 – 3 and 5 and maintaining pre-
scaled notation, we have

$\frac{\partial\boldsymbol{u}}{\partial
t}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}=-\nabla
p+\sqrt{\frac{Pr}{Ra}}\,\nabla^{2}\boldsymbol{u}-\frac{2}{Ro}\,\boldsymbol{\hat{z}}\times\boldsymbol{u}+\theta\,\boldsymbol{\hat{z}}+{\frac{1}{Ri_{b}}}\,\boldsymbol{\hat{x}},$
(10) $\nabla\cdot\boldsymbol{u}=0,$ (11) $\frac{\partial\theta}{\partial
t}+\boldsymbol{u}\cdot\nabla\theta=\sqrt{\frac{1}{Pr\cdot
Ra}}\,\nabla^{2}\theta\text{{, and}}$ (12) $v_{n}=\sqrt{\frac{1}{Pr\cdot
Ra\cdot St^{2}}}\,\boldsymbol{n}\cdot\nabla\theta\hskip
5.69046pt\text{at}\hskip 2.84544ptz=h(x,y,t),$ (13)

where the dimensionless temperature,

$\theta=\frac{T-T_{m}}{\Delta T},$ (14)

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, and are defined as

$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}{;}$
(15) $St=\frac{L_{s}}{C_{p}\,\Delta T}{;}\quad\mbox{and}\quad
Ro=\sqrt{\frac{g\,\alpha\,\Delta T}{\Omega^{2}\,H}},$ (16)

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 in 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

$Re=\sqrt{\frac{Ra}{Ri_{b}Pr}}.$ (17)

In cases where both rotation and shear are present, we define $\Sigma$ as the
ratio of centrifugal to applied body forces,

$\Sigma=\frac{Ri_{b}}{Ro^{2}}=\frac{\Omega^{2}H}{F_{p}},$ (18)

such that mean shear dominates when $\Sigma<1$.

The dimensionless temperature boundary conditions are:

$\theta(x,y,z=0,t)=1\quad\mbox{and}\quad\theta(x,y,z=H,t)=0,$ (19)

at the bottom and top surfaces, respectively, and

$\theta(x,y,z=h,t)=0,$ (20)

at the phase boundary. The no-slip and no-penetration boundary conditions for
the velocity field remain unaltered after the non-dimensionalisation.

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

$Nu={\sqrt{RaPrSt^{2}}\,h\,\frac{dh}{dt}},$ (21)

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 and Wettlaufer (2021),
although their Stefan number is the inverse of that defined here.

##  III Numerical method

Equations 11 – 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 extensively
validated for free-shear flows (Singhal  _et al._ , 2021a, b) and Rayleigh-
Bénard convection (Ravichandran  _et al._ , 2020; Ravichandran and Wettlaufer,
2020, 2021). 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
and Wettlaufer (2021). 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 equation 1 produces the known analytical plane-Poiseuille
velocity profile with an exponential relaxation timescale of
$(\pi^{2}\nu)^{-1}$ for the peak velocity.

##  IV Results

Here, we present our results on the effects of mean shear, buoyancy and
rotation on the evolution of the phase boundary. 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 and Wettlaufer,
2020). 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.

$Ro$ | $Pr$ | $Ri_{b}$ | $Re$ | $h_{0}$ | $(L_{x}\times L_{y})\times L_{z}$ | $N_{x},N_{y},N_{z}$ | $\Delta t$ | Figures
---|---|---|---|---|---|---|---|---
$\infty$ | $1$ | $\infty$ | $0$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  2, 7
$\infty$ | $1$ | $100$ | $35.4$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  2, 7
$\infty$ | $1$ | $10$ | $111.8$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  2, 3, 4, 6, 7
$\infty$ | $1$ | $4$ | $176.8$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  2, 3, 4, 6, 7
$\infty$ | $1$ | $1$ | $353.6$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 256$ | $2.8\times 10^{-4}$ |  2, 5, 6, 7
$Ro$ | $Pr$ | $Ri_{b}$ | $\Sigma$ | $h_{0}$ | $(L_{x}\times L_{y})\times L_{z}$ | $N_{x},N_{y},N_{z}$ | $\Delta t$ | Figures
$0.6325$ | $5$ | $100$ | $250$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  9, 15
$0.6325$ | $5$ | $10$ | $25$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  9, 15
$0.6325$ | $5$ | $1$ | $2.5$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  9, 10, 15
$0.6325$ | $5$ | $0.5$ | $1.25$ | $1$ | $16^{2}\times 2$ | $1024^{2}\times 128$ | $7\times 10^{-4}$ |  9, 15
$0.6325$ | $5$ | $0.2$ | $0.5$ | $0.5$ | $8^{2}\times 2$ | $512^{2}\times 128$ | $7\times 10^{-4}$ |  12, 13, 16
$0.6325$ | $5$ | $0.1$ | $0.25$ | $0.5$ | $8^{2}\times 2$ | $512^{2}\times 128$ | $7\times 10^{-4}$ |  11, 12, 13, 14, 16, 17
Table 1:  List of cases reported here. In all cases the Rayleigh number is
$Ra=1.25\times 10^{5}$, the Stefan number $St=1$ and the penalisation
parameter is $\eta=1.4\times 10^{-3}$. We also give the bulk Reynolds number
$Re$ for the non-rotating cases (equation 17) and the parameter $\Sigma$ for
the rotating cases (equation 18). The solutions obtained are independent of
the grid resolution, time step $\Delta t$ and volume penalisation parameter
$\eta$.

###  IV.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\left[1,\infty\right)$ in this
section, with the other dimensionless parameters fixed at $Ra=1.25\times
10^{5}$ and $Pr=St=1$.

In figures 2(a) – 2(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) and (d). The
mean liquid layer depths are comparable in (a-c), while the liquid depth in
(d) is about $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._ , 2019; Purseed  _et al._ , 2020) and three (Esfahani
_et al._ , 2018) 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) – 2(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 and Wettlaufer, 2019). Such features
have also been observed in the recent study of Couston  _et al._ (2020). 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, 2021).

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 the $Re$ is not too large (Clever and Busse,
1991). These longitudinal rolls can be discerned in the contours of the
vertical component 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.

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.

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._ (1980). This is due to the fact that in their experiments the $Re$ based
on the boundary-layer thickness was $O(10^{4})$, whereas in our simulations
the bulk $Re$ is $O(10^{2}-10^{3})$, and the wavelength of transverse patterns
is $O(10^{3})$ times the viscous sublayer thickness (Claudin  _et al._ , 2017)
and hence longer than our simulation domain. Longitudinal features are also
observed in flows over eroding (Allen, 1971) and dissolving boundaries (Cohen
_et al._ , 2016, 2020), where the associated grooves eventually merge to
produce scallops (see also Bushuk  _et al._ , 2019). 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, 2021).

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 a function of
$y$ and $t$ in figure 4. At $t=150$, for $Ri_{b}=10$ the six corrugations have
a wavelength of 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 4(a) and (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._
, 2018) and rotating cases (Ravichandran and Wettlaufer, 2021), where there
are there 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 $6$ domes of nearly
equal widths (of about $2.66$), and (b) $Ri_{b}=4$ showing the emergence of
$3$ 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) and (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 5(e) and (f)), which is homogeneous along the $x$ (figures
5(c). 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_
(c.f. 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}$. 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}=\left[\nu\left(\partial\bar{u}/\partial z\right)_{z=0}\right]^{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}$.

Figure 7:  (a) The average thickness of the solid layer $h_{s}$ (red curves,
left axis) and the melting Nusselt number (from eq. 21; blue curves, right
axis); and (b) the 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 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.

Figure 8:  The maximum values of (a) $Nu$ and (b) $\sigma_{h}$ in figure 7 as
functions of $Ri_{b}$.

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 the $Nu$
and 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 3-D flow. The effects are qualitatively consistent with those
observed in the study of pure Rayleigh-Bénard-Poiseuille flow in 3-D by
Scagliarini  _et al._ (2014).

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 homogenisation 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) and (b),
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 and Wettlaufer, 2021).

###  IV.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 §IV.1 and choose $Pr=5$ to give $Ro=0.6325$, which is fixed for all the
cases discussed in this section. The range of $Ri_{b}$ explored here is
$Ri_{b}\in\left[0.1,100\right]$, slightly larger than in §IV.1.

Figures 9(a)-(d) show the phase-boundary geometry 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, 1969; Ravichandran and Wettlaufer, 2021), 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 2(c), (d) and
9(c), (d)). The alignment of the longitudinal corrugations/grooves is 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 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$. Characteristic melting by rotating
convection (see e.g. Ravichandran and Wettlaufer, 2021) is seen in (a) 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 timescales 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 $\Sigma$ 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, figures 10(a,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 (figures
10(c-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 $\Sigma$ 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) and (f) to the melting morphology shown in figure 9(c)).

The loss of coherence of the columnar vortices is complete for $Ri_{b}\leq
0.2$, and the longitudinal rolls now clearly appear in the flow field. This is
seen in figure 11, which shows the contours of vertical velocity and
temperature for $Ri_{b}=0.1$.

Figure 11:  Contours of vertical velocity $w$ and temperature $\theta$ as in
figure 10, and for the same parameters except $Ri_{b}=0.1$ ($\Sigma=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.

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 the
one seen in thermal convection with mean shear flows, but without rotation
(Clever and Busse, 1991, 1992; Pabiou  _et al._ , 2005). These interfering
waves are clearly reflected in the patterns at the phase boundary, as seen in
figure 12(d).

Figure 12:  The 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, $\Sigma<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 $h_{0}=0.5$ in these simulations. The angle between the
grooves and the $x$-axis is smaller for smaller $Ri_{b}$.

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 the mean horizontal velocity makes with the
applied pressure gradient for $Ri_{b}=0.1$.

Figure 13:  (a) The horizontally averaged mean values of the velocities $u$
and $v$ along the $x$ and $y$ directions respectively, and (b) the 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), as
shown by the labels for $Ri_{b}$ in the inset. 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.

These plots show that despite the strength of the applied pressure gradient,
the flow in the bulk remains in geostrophic balance (and thus the angle made
by the flow relative to the $x$ direction is $\phi=-\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 $-\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 12(a) and (b)) and $Ri_{b}=0.1$ (figures 12(c) and
(d)). In contrast to the large $Ri_{b}$ cases, these grooves are neither
aligned parallel nor perpendicular to the direction of the applied pressure
gradient but at an intermediate angle.

The Ekman layer at the upper boundary may be defined as the region where the
angle $\phi\neq-\pi/2$ (see figure 13(b)). 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.

Figure 15:  (a) The thickness of the solid layer $h_{s}$ (red curves, left
axis) and the melting Nusselt number (from eq. 21; blue curves, right axis);
and (b) the roughness of the solid-liquid interface $\sigma_{h}$ for
$Ri_{b}=100,10,1,0.5$ (with $\Sigma=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$.

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 and Wettlaufer, 2021).

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.

Figure 16:  (a) The thickness of the solid layer and (b) the Nusselt number as
a function of time for $Ro=0.6325$, $Ra=1.25\times 10^{5}$, $Pr=5$, $St=1$,
with $\Sigma=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$.

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
_et al._ , 2015). 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$.

In geostrophic convection, the interaction between the thermal and momentum
boundary layers plays an important role (Rossby, 1969; Julien  _et al._ ,
2012; King  _et al._ , 2012). 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,
becoming 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).

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 boundary, while the dotted lines show the average height of
the fluid layer.

The different phase-boundary geometries suggest the existence of four distinct
regimes: (1) $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) $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 9(b,c) and 10(a,b)); (3)
$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)
$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 becoming sinusoidal.

##  V 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 (a) $Ro=\infty,1$; (b) $Ra=1.25\times 10^{5}$; (c)
$St=1$; (e) $Pr=1,5$; and (f) $Ri_{b}\in\left[0.1,\infty\right)$. 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.

  1. 1.

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 quasi-two-dimensional. This is seen
in figures 2, and 3(a,b).

  2. 2.

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 and Busse, 1991). This is in agreement with the numerical results of
Couston  _et al._ (2020), but in contrast to the experimental results of
Gilpin  _et al._ (1980). In our simulations, the bulk $Re=O(10^{2}-10^{3})$,
whereas in the experiments of Gilpin  _et al._ (1980) the $Re$ based on the
boundary-layer thickness is $O(10^{4})$, so that the flow regimes are not
directly comparable.

  3. 3.

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
direction of the applied pressure gradient. The dynamics is governed by the
balance between mean shear and rotation, and is represented by the parameter
$\Sigma$ (equation 18), with the imposed shear becoming dominant for
$\Sigma<1$.

  1. 1.

For $Ri_{b}=100$ ($\Sigma\gg 1$), the shear flow is weak and the phase-
boundary consists of dome-like features.

  2. 2.

On reducing $Ri_{b}$ to $10$ ($\Sigma=25$), these dome-like features begin to
merge. A further reduction in $Ri_{b}$ leads to the formation of grooves at
the phase boundary.

  3. 3.

For $Ri_{b}\in\left[0.2,0.5\right]$ ($\Sigma\in\left[0.5,1.25\right]$), the
merger is of the dome-like features is complete. The rolls develop a wave-like
instability for $Ri_{b}=0.2$ and $0.1$ ($\Sigma=0.5$ and $0.25$,
respectively), which leads to sinusoidal evolution of the phase boundary. This
instability is qualitatively similar to that observed in thermal convection in
the presence of a mean shear flow (Clever and Busse, 1991, 1992; Pabiou  _et
al._ , 2005).

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 and Jamtveit, 2009; Neufeld  _et al._ ,
2010; Bushuk  _et al._ , 2019; Weady  _et al._ , 2022, and refs. 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._ , 2012;
Julien  _et al._ , 2012). 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.

## Acknowledgements

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.

## References

  * Epstein and Cheung (1983) M. Epstein and F. B. Cheung, Ann. Rev. Fl. Mech. 15, 293 (1983).
  * Glicksman  _et al._ (1986) M. E. Glicksman, S. R. Coriell, and G. B. McFadden, Annu. Rev. Fl. Mech. 18, 307 (1986).
  * Huppert (1986) H. E. Huppert, J. Fluid Mech. 173, 557 (1986).
  * Worster (2000) M. G. Worster, in  _Perspectives in Fluid Dynamics — a Collective Introduction to Current Research_ , edited by G. Batchelor, H. Moffatt, and M. Worster (Cambridge University Press, 2000) pp. 393 – 446\.
  * Meakin and Jamtveit (2009) P. Meakin and B. Jamtveit, Proc. R. Soc. A 466, 659 (2009).
  * Hewitt (2020) I. J. Hewitt, Annu. Rev. Fl. Mech. 52, 145 (2020).
  * Davis  _et al._ (1984) S. H. Davis, U. Müller, and C. Dietsche, J. Fluid Mech. 144, 133 (1984).
  * Dietsche and Müller (1985) C. Dietsche and U. Müller, J. Fluid Mech. 161, 249 (1985).
  * Wettlaufer  _et al._ (1997) J. S. Wettlaufer, M. G. Worster, and H. E. Huppert, J. Fluid Mech. 344, 291 (1997).
  * Worster (1997) M. G. Worster, Ann. Rev. Fl. Mech. 29, 91 (1997).
  * Davies Wykes  _et al._ (2018) M. S. Davies Wykes, J. M. Huang, G. A. Hajjar, and L. Ristroph, Phys. Rev. Fluids 3, 043801 (2018).
  * Delves (1968) R. T. Delves, J. Cryst. Growth 3, 562 (1968).
  * Delves (1971) R. T. Delves, J. Cryst. Growth 8, 13 (1971).
  * Gilpin  _et al._ (1980) R. R. Gilpin, T. Hirata, and K. C. Cheng, J. Fluid Mech. 99, 619 (1980).
  * Coriell  _et al._ (1984) S. R. Coriell, G. B. McFadden, R. F. Boisvert, and R. F. Sekerka, J. Cryst. Growth 69, 15 (1984).
  * Forth and Wheeler (1989) S. A. Forth and A. A. Wheeler, J. Fluid Mech. 202, 339 (1989).
  * Feltham and Worster (1999) D. L. Feltham and M. G. Worster, J. Fluid Mech. 391, 337 (1999).
  * Neufeld and Wettlaufer (2008a) J. A. Neufeld and J. S. Wettlaufer, J. Fluid Mech. 612, 363 (2008a).
  * Neufeld and Wettlaufer (2008b) J. A. Neufeld and J. S. Wettlaufer, J. Fluid Mech. 612, 339 (2008b).
  * Dallaston  _et al._ (2015) M. C. Dallaston, I. J. Hewitt, and A. J. Wells, J. Fluid Mech. 785, 109 (2015).
  * Ramudu  _et al._ (2016) E. Ramudu, B. H. Hirsh, P. Olson, and A. Gnanadesikan, J. Fluid Mech. 798, 572 (2016).
  * Bushuk  _et al._ (2019) M. Bushuk, D. M. Holland, T. P. Stanton, A. Stern, and C. Gray, J. Fluid Mech. 873, 942 (2019).
  * McPhee (2008) M. G. McPhee,  _Air-Ice-Ocean Interaction: Turbulent Ocean Boundary Layer Exchange Processes_ (Springer, New York, NY, 2008).
  * Kurz and Fischer (1984) W. Kurz and D. J. Fischer,  _Fundamentals of Solidification_ (Trans Tech Publications, Switzerland-Germany-UK- USA, 1984).
  * Toppaladoddi and Wettlaufer (2019) S. Toppaladoddi and J. S. Wettlaufer, J. Fluid Mech. 868, 648 (2019).
  * Worster  _et al._ (2021) M. G. Worster, S. S. L. Peppin, and J. S. Wettlaufer, J. Fluid Mech. 914, A28 (2021).
  * Worster (1992) M. G. Worster, J. Fluid Mech. 237, 649 (1992).
  * Drazin and Reid (2004) P. G. Drazin and W. H. Reid,  _Hydrodynamic Stability_ (Cambridge University Press, Cambridge, 2004).
  * Neufeld  _et al._ (2006) J. A. Neufeld, J. S. Wettlaufer, D. L. Feltham, and M. G. Worster, J. Fluid Mech. 549, 442 (2006).
  * Sample and Hellawell (1982) A. Sample and A. K. Hellawell, Metall. Trans. B 13, 495 (1982).
  * Sample and Hellawell (1984) A. K. Sample and A. Hellawell, Metall. Trans. A 15, 2163 (1984).
  * Lu and Chen (1997) J. W. Lu and F. Chen, Int. J. Heat Mass Transfer 40, 237 (1997).
  * Chung and Chen (2000) C. A. Chung and F. Chen, J. Fluid Mech. 412, 93 (2000).
  * Chung and Chen (2003) C. A. Chung and F. Chen, J. Fluid Mech. 477, 381 (2003).
  * Esfahani  _et al._ (2018) B. R. Esfahani, S. C. Hirata, S. Berti, and E. Calzavarini, Phys. Rev. Fluids 3, 053501 (2018).
  * Favier  _et al._ (2019) B. Favier, J. Purseed, and L. Duchemin, J. Fluid Mech. 858, 437 (2019).
  * Purseed  _et al._ (2020) J. Purseed, B. Favier, L. Duchemin, and E. W. Hester, Phys. Rev. Fluids 5, 023501 (2020).
  * Toppaladoddi (2021) S. Toppaladoddi, J. Fluid Mech. 919, A28 (2021).
  * Hirata  _et al._ (1979a) T. Hirata, R. R. Gilpin, and K. C. Cheng, Int. J. Heat Mass Transfer 22, 1435 (1979a).
  * Hirata  _et al._ (1979b) T. Hirata, R. R. Gilpin, K. C. Cheng, and E. M. Gates, Int. J. Heat Mass Transfer 22, 1425 (1979b).
  * Couston  _et al._ (2020) L.-A. Couston, E. Hester, B. Favier, J. R. Taylor, P. R. Holland, and A. Jenkins, J. Fluid Mech. 911, A44 (2020).
  * Ravichandran and Wettlaufer (2021) S. Ravichandran and J. Wettlaufer, J. Fluid Mech. 916, A28 (2021).
  * Rossby (1969) H. T. Rossby, J. Fluid Mech. 36, 309 (1969).
  * Wells and Wettlaufer (2008) M. G. Wells and J. S. Wettlaufer, Geophys. Res. Lett. 35, L03501 (2008).
  * Singhal  _et al._ (2021a) R. Singhal, S. Ravichandran, S. S. Diwan, and G. L. Brown, in  _Proceedings of 16th Asian Congress of Fluid Mechanics_ (Springer, 2021) pp. 403–411.
  * Singhal  _et al._ (2021b) R. Singhal, S. Ravichandran, R. Govindarajan, and S. S. Diwan, arXiv:2103.16415 (2021b).
  * Ravichandran  _et al._ (2020) S. Ravichandran, E. Meiburg, and R. Govindarajan, J. Fluid Mech. 899, A27 (2020).
  * Ravichandran and Wettlaufer (2020) S. Ravichandran and J. S. Wettlaufer, J. Fluid Mech. 897, A24 (2020).
  * Clever and Busse (1991) R. M. Clever and F. H. Busse, J. Fluid Mech. 229, 517 (1991).
  * Claudin  _et al._ (2017) P. Claudin, O. Durán, and B. Andreotti, J. Fluid Mech. 832, R2 (2017).
  * Allen (1971) J. R. L. Allen, J. Fluid Mech. 49, 49 (1971).
  * Cohen  _et al._ (2016) C. Cohen, M. Berhanu, J. Derr, and S. C. D. Pont, Phys. Rev. Fluids 1, 50508 (2016).
  * Cohen  _et al._ (2020) C. Cohen, M. Berhanu, J. Derr, and S. C. D. Pont, Phys. Rev. Fluids 5, 053802 (2020).
  * Scagliarini  _et al._ (2014) A. Scagliarini, Á. Gylfason, and F. Toschi, Phys. Rev. E 89, 043012 (2014).
  * Clever and Busse (1992) R. M. Clever and F. H. Busse, J. Fluid Mech. 234, 511 (1992).
  * Pabiou  _et al._ (2005) H. Pabiou, S. Mergui, and C. Bénard, J. Fluid Mech. 542, 175 (2005).
  * Toppaladoddi  _et al._ (2015) S. Toppaladoddi, S. Succi, and J. S. Wettlaufer, EPL 111, 44005 (2015).
  * Julien  _et al._ (2012) K. Julien, E. Knobloch, A. M. Rubio, and G. M. Vasil, Phys. Rev. Lett. 109, 254503 (2012).
  * King  _et al._ (2012) E. M. King, S. Stellmach, and J. M. Aurnou, J. Fluid Mech. 691, 568 (2012).
  * Neufeld  _et al._ (2010) J. A. Neufeld, R. E. Goldstein, and M. G. Worster, J. Fluid Mech. 647, 287 (2010).
  * Weady  _et al._ (2022) S. Weady, J. Tong, A. Zidovska, and L. Ristroph, Phys. Rev. Lett. 128, 044502 (2022).