1 Introduction
Flow of a melt over its solid phase can profoundly influence the latter’s evolution and stability (e.g. Epstein & Cheung Reference Epstein and Cheung1983; Glicksman, Coriell & McFadden Reference Glicksman, Coriell and McFadden1986; Davis Reference Davis1990; Huppert Reference Huppert1990; Worster Reference Worster, Batchelor, Moffatt and Worster2000). Examples abound in both natural (e.g. Untersteiner & Badgley Reference Untersteiner and Badgley1965; Wettlaufer Reference Wettlaufer1991; McPhee Reference McPhee2008; Meakin & Jamtveit Reference Meakin and Jamtveit2009; Solari & Parker Reference Solari and Parker2013; Ramudu et al. Reference Ramudu, Hirsh, Olson and Gnanadesikan2016; Claudin, Durán & Andreotti Reference Claudin, Durán and Andreotti2017) and engineering (e.g. Delves Reference Delves1968, Reference Delves1971; Forth & Wheeler Reference Forth and Wheeler1989) settings. Flows over phase-changing boundaries can be grouped into the following two categories: (i) free flows, which arise due to density differences created during solidification (Davis, Müller & Dietsche Reference Davis, Müller and Dietsche1984; Liu, Ning & Ecke Reference Liu, Ning and Ecke1996; Wettlaufer, Worster & Huppert Reference Wettlaufer, Worster and Huppert1997; Worster Reference Worster1997; Davies Wykes et al. Reference Davies Wykes, Huang, Hajjar and Ristroph2018), and (ii) forced flows, which are typically shear driven, and are introduced to control morphological and/or hydrodynamical instabilities (Delves Reference Delves1968, Reference Delves1971; Coriell et al. Reference Coriell, McFadden, Boisvert and Sekerka1984; Forth & Wheeler Reference Forth and Wheeler1989).
In the absence of an external flow, the rates of freezing are typically sufficiently large so that a planar solid–binary liquid interface will become highly convoluted, leading to one of the two components being trapped in the interstices of the crystals of the other component (Worster Reference Worster, Batchelor, Moffatt and Worster2000). In engineering, the imposition of a flow was motivated by controlling the instability, whereas, in natural settings, it is often an unavoidable part of the environment (e.g. Delves Reference Delves1968, Reference Delves1971; Coriell et al. Reference Coriell, McFadden, Boisvert and Sekerka1984; Forth & Wheeler Reference Forth and Wheeler1989; Schulze & Davis Reference Schulze and Davis1994, Reference Schulze and Davis1995, Reference Schulze and Davis1996; Feltham & Worster Reference Feltham and Worster1999; Feltham, Worster & Wettlaufer Reference Feltham, Worster and Wettlaufer2002; Neufeld et al. Reference Neufeld, Wettlaufer, Feltham and Worster2006; Neufeld & Wettlaufer Reference Neufeld and Wettlaufer2008a ,Reference Neufeld and Wettlaufer b ; Camporeale & Ridolfi Reference Camporeale and Ridolfi2012). Here, we focus on understanding the effects of shear and buoyancy on the directional solidification of a pure melt, but we shall review the results on binary mixtures as well, because there are some commonalities in the dynamics of the two systems. The principal difference between the two systems is that in a binary mixture the evolving morphology affects the growth of disturbances at the mush–liquid interface through fluid flows in both the mushy layer interior and the liquid phase. However, no such mechanism exists in the pure case.
Some of the first studies to investigate the effects of shear-driven flows on directional solidification of a binary alloy using linear stability analysis are those of Delves (Reference Delves1968, Reference Delves1971) and Coriell et al. (Reference Coriell, McFadden, Boisvert and Sekerka1984). Delves (Reference Delves1968, Reference Delves1971) studied the effects of a parabolic flow on morphological instability and found that the flow suppresses the instability, with the degree of suppression depending on the material considered. He also found that the flow gives rise to travelling waves along the interface. Coriell et al. (Reference Coriell, McFadden, Boisvert and Sekerka1984) studied the effects of Couette flow on the morphological and thermosolutal instabilities during directional solidification of a lead–tin alloy. Their findings suggest that Couette flow suppresses the onset of thermosolutal instability to a larger degree than the onset of morphological instability. However, the use of Couette flow as the base-state velocity profile seems incompatible with the momentum-balance equations (Coriell et al. Reference Coriell, Cordes, Boettinger and Sekerka1980), which admit the asymptotic-suction boundary-layer profile (Drazin & Reid Reference Drazin and Reid2004) as their solution.
Forth & Wheeler (Reference Forth and Wheeler1989) studied directional solidification of a binary alloy in the presence of an asymptotic-suction boundary-layer flow. They focused on (i) understanding how the fluid flow affects the morphological instability, and (ii) understanding how the freezing interface affects the shear-flow instability. Under certain conditions they find that the shear flow only leads to the generation of travelling waves along the interface, and that the speed of these waves varies linearly with the imposed flow speed. However, under the same conditions, the freezing interface was found to have negligible effects on the hydrodynamic instability.
The structure resulting from the instability of the solid–binary liquid interface is known as a mushy layer (Worster Reference Worster, Batchelor, Moffatt and Worster2000), and is modelled as a chemically reacting porous medium (Worster Reference Worster1991). The most common example of mushy layer is the sea ice found in Earth’s polar regions (Feltham et al. Reference Feltham, Untersteiner, Wettlaufer and Worster2006). Here, compositional convection can be induced both in the mushy layer, which contains brine trapped between ice crystals, and in the sea water, which is gravitationally unstable due to high concentration of salt – rejected during solidification – close to the ice–water interface (Worster Reference Worster1992). These modes of convection are termed mushy and boundary-layer modes, respectively (Worster Reference Worster1992), and have been observed in the laboratory (Wettlaufer et al. Reference Wettlaufer, Worster and Huppert1997). It is intuitive that in the presence of a shear flow, the evolution of any incipient perturbation at the mush–liquid interface should depend on the interaction between the flows in the melt and mushy layer.
By neglecting the effects of buoyancy in both the bulk melt and the mushy layer, Feltham & Worster (Reference Feltham and Worster1999) investigated the effects of forced flow of inviscid and viscous melts on the morphology of a mushy layer. They found that an external flow over a corrugated mush–liquid interface results in a pressure perturbation along the interface that drives flow in the mushy layer, and under certain conditions this leads to the growth of the perturbations with a wavelength commensurate with the depth of the mushy layer. The perturbed heat flux from the liquid was found to have no influence on the evolution of the perturbation and was only responsible for introducing travelling waves at the interface.
Neufeld & Wettlaufer (Reference Neufeld and Wettlaufer2008a ,Reference Neufeld and Wettlaufer b ) studied the effects of shear flow on the mushy- and boundary-layer modes of convection using both theory and experiments. They found that; (i) below a critical value of the shear-flow velocity, both modes of convection are moderately suppressed; (ii) above a critical shear-flow velocity, the stability of both modes of convection decreases monotonically with the strength of the flow; (iii) for sufficiently strong shear flow, striations of zero solid fraction transverse to the flow direction are generated. These striations are quasi-two-dimensional and form because of localized dissolution and growth of the mushy layer, which in turn are due to the interplay between shear and buoyancy.
Relative to binary mixtures, there have been far fewer studies of the influence of external flows on the directional solidification of pure melts. One of the first experimental studies was by Gilpin, Hirata & Cheng (Reference Gilpin, Hirata and Cheng1980), who investigated the evolution of a layer of pure ice in contact with a turbulent flow in a closed-loop water tunnel with an upper free surface. A layer of ice was formed on a surface that was held at a temperature less than the melting temperature and a shear flow was maintained over the ice layer, with the far-field temperature greater than the melting point. Before starting the flow, the ice–water interface was perturbed by melting a groove into the ice layer. Under certain conditions, the perturbation at the ice–water interface was observed to grow, leading to the formation of a ‘rippled’ surface. They found that the heat transfer rate over the rippled surface was
$30$
–
$60\,\%$
larger than that on a planar surface and the evolution of the ice layer was wholly attributed to the overlying shear flow. However, one crucial point that Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) evidently overlooked is that because the far-field temperature of water was greater than the melting point, the water column above the ice layer was unstably stratified due to the
$4\,^{\circ }\text{C}$
density maximum, which can exert a controlling influence on heat flux (Veronis Reference Veronis1963; Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2018). Indeed, an estimate of the Rayleigh number – defined as the ratio of buoyancy to viscous forces – for the experiments of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) gives the value as
${\approx}8\times 10^{7}$
, which is indicative of the important role of buoyancy.
Here, motivated in part by the experiments of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980), we study the effects of shear and buoyancy on the phase evolution of a pure melt. Specifically, we study the growth of perturbations at a solid–liquid interface in the presence of Couette flow and Rayleigh–Bénard convection using linear stability analysis. The reason for our choice of the Rayleigh–Bénard–Couette system is to have an analytically tractable system where the relative effects of shear and buoyancy on the stability of the phase boundary can be discerned more easily than in a more complex flow.
2 Governing equations
To perform a linear stability analysis, we consider the domain shown in figure 1. The length of the cell is
$L_{x}$
and the depth of the cell is
$L_{z}$
. At the initial instant the solid occupies the region
$h_{0}\leqslant z\leqslant L_{z}$
, and the liquid occupies
$0\leqslant z\leqslant h_{0}$
. The solid–liquid interface is planar and is at
$z=h_{0}$
. The initial thickness of the solid layer is
$d_{0}$
, and hence
$L_{z}=h_{0}+d_{0}$
. The upper surface is maintained at a temperature
$T_{c}$
and the lower surface is maintained at
$T_{h}$
. The temperatures are such that
$T_{h}>T_{m}>T_{c}$
, where
$T_{m}$
is the melting temperature of the solid. The liquid considered has a linear equation of state, hence the liquid column in unstably stratified. The bottom surface moves at a constant horizontal velocity
$U_{\infty }$
, as shown in figure 1.
Schematic of the domain considered here.

The governing equations in the different regions are as follows.
2.1 Liquid
The continuity, Boussinesq, and heat-balance equations are
respectively. Here,
$\boldsymbol{u}(\boldsymbol{x},t)=(u,v,w)$
is the velocity field,
$\unicode[STIX]{x1D70C}_{0}$
is the reference density,
$p(\boldsymbol{x},t)$
is the pressure field,
$g$
is acceleration due to gravity,
$\unicode[STIX]{x1D6FC}$
is the thermal expansion coefficient,
$T_{l}(\boldsymbol{x},t)$
is the temperature field,
$\unicode[STIX]{x1D708}$
is the kinematic viscosity and
$\unicode[STIX]{x1D705}$
is the thermal diffusivity. To simplify matters, we assume the liquid and solid phases have the same density (
$\unicode[STIX]{x1D70C}_{0}$
) and thermal diffusivity (
$\unicode[STIX]{x1D705}$
).
2.2 Solid
The temperature field in the solid,
$T_{s}(\boldsymbol{x},t)$
, is governed by diffusion, viz.
2.3 Solid–liquid interface
At the solid–liquid interface, we have the Stefan condition
where
$L_{s}$
is the latent heat of fusion,
$\boldsymbol{n}$
is the unit vector pointing into the liquid,
$\boldsymbol{q}_{\boldsymbol{s}}=-k\,\unicode[STIX]{x1D735}T_{s}|_{z=h^{+}}$
is the heat flux away from the interface into the solid and
$\boldsymbol{q}_{\boldsymbol{l}}=-k\,\unicode[STIX]{x1D735}T_{l}|_{z=h^{-}}$
is the heat flux towards the interface from the liquid.
2.4 Boundary conditions
The boundary conditions for heat equation in the solid are
and those for the advection–diffusion equation in the liquid are
The velocity field satisfies
and
We non-dimensionalize these equations by choosing
$U_{\infty }$
as the velocity scale;
$h_{0}$
as the length scale,
$t_{0}=h_{0}^{2}/\unicode[STIX]{x1D705}$
as the time scale,
$p_{0}=\unicode[STIX]{x1D70C}_{0}\,U_{\infty }\,\unicode[STIX]{x1D705}/h_{0}$
as the pressure scale and
$\unicode[STIX]{x0394}T=T_{h}-T_{m}$
as the temperature scale. Using these in (2.2)–(2.5), and maintaining the pre-scaled notation, we have
and
where,
There are five governing parameters, which are
where,
$Ra$
,
$Pe$
,
$Pr$
and
${\mathcal{S}}$
are the Rayleigh, Péclet, Prandtl and Stefan numbers, respectively. The ratio of the temperature differences across the liquid and the solid regions is denoted by
$\unicode[STIX]{x1D6EC}$
.
The thermal and velocity boundary conditions now become
3 Linear stability analysis
We now perform linear stability analysis on (2.10)–(2.14), with the boundary conditions (2.18)–(2.21).
3.1 Base-state solutions
All variables in the base state are assumed to be steady and horizontally homogeneous.
3.1.1 Liquid
The base-state velocity and temperature profiles are taken to be
$u^{(0)}(z)$
and
$\unicode[STIX]{x1D703}_{l}^{(0)}(z)$
. Solving the equations of motion subject to the boundary conditions gives
3.1.2 Solid
The solution to the heat equation for the base-state temperature field in the solid is given by
3.1.3 Interface
In the base state, we assume that the heat fluxes away from and towards the interface balance, so that the initial thickness of the solid layer is constant. Hence, the Stefan condition is
which gives
3.2 Equations for the perturbation amplitudes
We introduce a normal mode perturbation of the interface given by
This in turn leads to perturbations in the liquid and solid layers so that the total velocity, pressure and temperature fields become
$$\begin{eqnarray}\left[\begin{array}{@{}c@{}}u(x,y,z,t)\\ v(x,y,z,t)\\ w(x,y,z,t)\\ p(x,y,z,t)\\ \unicode[STIX]{x1D703}_{l}(x,y,z,t)\\ \unicode[STIX]{x1D703}_{s}(x,y,z,t)\end{array}\right]=\left[\begin{array}{@{}c@{}}u^{(0)}(z)\\ 0\\ 0\\ p^{(0)}(z)\\ \unicode[STIX]{x1D703}_{l}^{(0)}(z)\\ \unicode[STIX]{x1D703}_{s}^{(0)}(z)\end{array}\right]+\unicode[STIX]{x1D716}\left[\begin{array}{@{}c@{}}\widehat{u}(z)\\ \widehat{v}(z)\\ \widehat{w}(z)\\ \widehat{p}(z)\\ \widehat{\unicode[STIX]{x1D703}_{l}}(z)\\ \widehat{\unicode[STIX]{x1D703}_{s}}(z)\end{array}\right]\exp (\text{i}kx+\text{i}my+\unicode[STIX]{x1D70E}t_{1}),\end{eqnarray}$$
where
$t_{1}=t/{\mathcal{S}}$
. The range of
${\mathcal{S}}$
in the experiments of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) was
${\mathcal{S}}\approx [4,7]$
, and hence we are interested in the limit
${\mathcal{S}}\gg 1$
, in which case the rate-controlling process is the release of latent heat, wherein the dynamics in the solid and liquid regions becomes quasi-steady (e.g. Feltham & Worster Reference Feltham and Worster1999).
3.2.1 Liquid
Linearizing equations (2.11)–(2.12) and using (3.6), we obtain the following equations for the amplitudes:
where
$D\equiv \text{d}/\text{d}z$
and
$\unicode[STIX]{x1D6FE}^{2}=k^{2}+m^{2}$
. The boundary conditions become
and
We now obtain a single equation for
$\widehat{\unicode[STIX]{x1D703}_{l}}$
. Following Forth & Wheeler (Reference Forth and Wheeler1989), we eliminate
$\widehat{u}$
and
$\widehat{v}$
from (3.8) and (3.9) to obtain
Eliminating
$\widehat{p}$
from (3.14) and (3.10) we obtain
Finally, eliminating
$\widehat{w}$
from (3.11) and (3.15) gives the following sixth-order ordinary differential equation for
$\widehat{\unicode[STIX]{x1D703}_{l}}$
$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle PrD^{6}\widehat{\unicode[STIX]{x1D703}_{l}}-[3Pr\unicode[STIX]{x1D6FE}^{2}+\text{i}ku^{(0)}Pe(1+Pr)]D^{4}\widehat{\unicode[STIX]{x1D703}_{l}}+(4\text{i}kPrPe)D^{3}\widehat{\unicode[STIX]{x1D703}_{l}}\nonumber\\ \displaystyle & & \displaystyle +\,[3Pr\unicode[STIX]{x1D6FE}^{4}+2\text{i}k\unicode[STIX]{x1D6FE}^{2}u^{(0)}Pe(1+Pr)-k^{2}Pe^{2}(u^{(0)})^{2}]D^{2}\widehat{\unicode[STIX]{x1D703}_{l}}\nonumber\\ \displaystyle & & \displaystyle -\,(4\text{i}k\unicode[STIX]{x1D6FE}^{2}PrPe-2k^{2}Pe^{2}u^{(0)})D\widehat{\unicode[STIX]{x1D703}_{l}}\nonumber\\ \displaystyle & & \displaystyle -\,[Pr\unicode[STIX]{x1D6FE}^{6}+\text{i}k\unicode[STIX]{x1D6FE}^{4}Peu^{(0)}(1+Pr)-k^{2}Pe^{2}\unicode[STIX]{x1D6FE}^{2}(u^{(0)})^{2}-RaPr\unicode[STIX]{x1D6FE}^{2}]\widehat{\unicode[STIX]{x1D703}_{l}}.\end{eqnarray}$$
The boundary conditions at
$z=0$
are
and
and those at
$z=1$
are
and
Equation (3.16), along with boundary conditions (3.17) and (3.18), is solved numerically using Chebfun (Driscoll, Bornemann & Trefethen Reference Driscoll, Bornemann and Trefethen2008).
3.2.2 Solid
The equation for
$\unicode[STIX]{x1D703}_{s}^{\prime }$
is
Using normal modes
$\unicode[STIX]{x1D703}_{s}^{\prime }=\widehat{\unicode[STIX]{x1D703}_{s}}\exp (\text{i}kx+\text{i}my+\unicode[STIX]{x1D70E}t_{1})$
, we have
with
as the boundary conditions. Equation (3.20) has solution
where
3.2.3 Interface
At
$O(\unicode[STIX]{x1D716})$
the Stefan condition becomes
from which it is evident that the heat flux from the liquid has considerable influence on the stability of the interface. Thus, generation of fluid motions with appreciable vertical velocities can lead to a larger perturbed heat flux, thereby making the interface unstable.
4 Results and discussion
4.1 Phase change with no shear flow
On setting
$Pe=0$
the present problem reduces to that of phase change in the presence of an unstably stratified column of liquid, which has been studied by Davis et al. (Reference Davis, Müller and Dietsche1984). When
$Pe=0$
, equation (3.16) is independent of
$Pr$
and hence so too is the critical Rayleigh number,
$Ra_{c}$
, at which convective motions develop (Davis et al.
Reference Davis, Müller and Dietsche1984; Chandrasekhar Reference Chandrasekhar2013). However, as shown by Davis et al. (Reference Davis, Müller and Dietsche1984),
$Ra_{c}$
and the critical wavenumber,
$\unicode[STIX]{x1D6FE}_{c}$
, are functions of
$d_{0}$
.
In figures 2 and 3 we compare
$Ra_{c}$
and
$\unicode[STIX]{x1D6FE}_{c}$
as functions of
$d_{0}$
with the calculations of Davis et al. (Reference Davis, Müller and Dietsche1984). The decrease in
$Ra_{c}$
with increasing
$d_{0}$
is due to the fact that the velocity boundary condition at the top surface for the liquid is ‘relaxed’ due to the presence of the moving boundary. The calculations of Davis et al. (Reference Davis, Müller and Dietsche1984) were focused on experiments using cyclohexane, for which we estimate
${\mathcal{S}}\approx 6-8$
, showing good agreement with our calculations in figures 2 and 3.
Comparison of
$Ra_{c}(d_{0})$
with Davis et al. (Reference Davis, Müller and Dietsche1984). Circles are the values from the present calculations, and the solid line is from Davis et al. (Reference Davis, Müller and Dietsche1984).

Comparison of
$\unicode[STIX]{x1D6FE}_{c}(d_{0})$
with Davis et al. (Reference Davis, Müller and Dietsche1984). Circles are the values from the present calculations, and the solid line is from Davis et al. (Reference Davis, Müller and Dietsche1984).

4.2 Effects of shear
4.2.1 Roll structure and its dependence on shear and the perturbation wave vector
The effects of shear flow on the perturbations at the interface depend on how the flow is aligned with respect to the perturbation wave vector
$\boldsymbol{\unicode[STIX]{x1D6FE}}=(k,m)$
(Chung & Chen Reference Chung and Chen2001; Neufeld & Wettlaufer Reference Neufeld and Wettlaufer2008b
). This dependence can be understood by following Chung & Chen (Reference Chung and Chen2001) and performing a Squire transformation of the base-state velocity. In our notation this is
Figures 4–6 show the perturbed temperature field for
$Ra=1700$
,
$Pe=0.5$
,
$\unicode[STIX]{x1D6FE}=3.021$
, and
$k=\unicode[STIX]{x1D6FE},m=0$
(figure 4),
$k=m=\unicode[STIX]{x1D6FE}/\sqrt{2}$
(figure 5) and
$k=0,m=\unicode[STIX]{x1D6FE}$
(figure 6), respectively. It is clearly seen that when
$u_{sq}^{(0)}=u^{(0)}(m=0)$
the rolls are aligned such that their axes are perpendicular to the direction of the flow; when
$0<u_{sq}^{(0)}<u^{(0)}(m\neq 0)$
the roll axes are aligned at a certain angle with the shear flow; and when
$u_{sq}^{(0)}=0$
$(m=\unicode[STIX]{x1D6FE})$
the roll axes are parallel to the shear flow.
A closer examination of figure 6 reveals that when
$k=0$
, the roll structure is completely unaffected by shear. Hence, shear has no effect on perturbations with wave vectors perpendicular to it (Chung & Chen Reference Chung and Chen2001; Neufeld & Wettlaufer Reference Neufeld and Wettlaufer2008b
). Noting this dependence on
$\boldsymbol{\unicode[STIX]{x1D6FE}}=(k,m)$
, we discuss the results in terms of
$\unicode[STIX]{x1D6FE}$
.
Perturbed temperature field for
$Ra=1700$
,
$Pe=0.5$
,
$\unicode[STIX]{x1D6FE}=3$
and
$k=\unicode[STIX]{x1D6FE},m=0$
in (a)
$x{-}z$
plane and (b)
$y{-}z$
plane. This case corresponds to
$u_{sq}^{(0)}=u^{(0)}$
. The dashed line denotes the solid–liquid interface.

Perturbed temperature field for
$Ra=1700$
,
$Pe=0.5$
,
$\unicode[STIX]{x1D6FE}=3$
and
$k=m=\unicode[STIX]{x1D6FE}/\sqrt{2}$
in (a)
$x{-}z$
plane and (b)
$y{-}z$
plane. This case corresponds to
$u_{sq}^{(0)}<u^{(0)}$
. The dashed line denotes the solid–liquid interface.

Perturbed temperature field for
$Ra=1700$
,
$Pe=0.5$
,
$\unicode[STIX]{x1D6FE}=3$
and
$k=0,m=\unicode[STIX]{x1D6FE}$
in (a)
$x{-}z$
plane and (b)
$y{-}z$
plane. This case corresponds to
$u_{sq}^{(0)}=0$
. The dashed line denotes the solid–liquid interface.

4.2.2 Effects on the instability
To understand the effects of shear on the instability of the convective flow, we solve equation (3.16) with
$Pe=0,0.15,0.5,2$
and
$5$
and a supercritical
$Ra$
of
$1700$
. Figure 7 shows the dispersion curve for the real part of the growth rate (
$\unicode[STIX]{x1D70E}_{r}$
) for
$d_{0}=0.1$
and different
$Pe$
. In the absence of shear, the most unstable mode has
$\unicode[STIX]{x1D6FE}=3.021$
, and clearly the interfacial instability is suppressed as the strength of the shear flow increases, with all modes decaying when
$Pe$
is as small as
$0.5$
. For example, in figure 8, we see that the growth rate becomes negative for
$Pe\approx 0.22$
, and asymptotes for
$Pe\geqslant 1$
.
Real growth rates
$\unicode[STIX]{x1D70E}_{r}$
as a function of wavenumber
$\unicode[STIX]{x1D6FE}$
when
$d_{0}=0.1$
and
$Ra=1700$
for different
$Pe$
. Shear has a strong stabilizing effect on the instability of the phase boundary.

The dependence of
$\unicode[STIX]{x1D70E}_{r}$
on
$Pe$
for
$d_{0}=0.1$
,
$Ra=1700$
and
$\unicode[STIX]{x1D6FE}=3.021$
, the most rapidly growing mode in the absence of shear flow (
$Pe=0$
). The growth rate becomes negative for
$Pe\geqslant 0.22$
. (See dotted red line.)

The introduction of shear flow leads to the stabilization of the interface, which is evidenced by the smaller values of
$\unicode[STIX]{x1D70E}_{r}$
relative to those for purely convective flow, and when
$\unicode[STIX]{x1D70E}_{r}>0$
we find travelling waves along the solid–liquid interface in the direction of the shear flow. As shown in figure 9, the
$\unicode[STIX]{x1D70E}_{i}(\unicode[STIX]{x1D6FE},Pe)$
curves display non-monotonic behaviour. This is because in the absence of shear flow, there are no travelling waves and the convective rolls are undistorted. Thus, for small
$Pe$
, these rolls are advected by the shear flow with little or no distortion. However, as
$Pe$
increases the convective and shear motions interact, leading to the excitation of a larger set of wavenumbers. This causes the convective rolls to lose their structural coherence.
Dispersion curves for
$\unicode[STIX]{x1D70E}_{i}$
with
$d_{0}=0.1,Ra=1700$
and different
$Pe$
.

These shear effects can be seen in figure 10, which shows the perturbed temperature field in the liquid and solid regions as a function of
$Pe$
.
We should note here that the values of
$\unicode[STIX]{x1D70E}$
in the experiments of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) may have an additional spatial dependence: because the flow is composed of both shear- and buoyancy-driven components, the turbulent flow field is spatially inhomogeneous. Hence, a perturbation originating at a particular location at the interface may have the magnitude and/or sign of its growth rate modified as it propagates along the interface. However, a theoretical study of the linear stability of the system avoids this complication.
Perturbed temperature field in the liquid and solid regions for
$\unicode[STIX]{x1D6FE}=3.021~(k=m=\unicode[STIX]{x1D6FE}/\sqrt{2})$
when
$d_{0}=0.1$
and
$Ra=1700$
. (a)
$Pe=0.0$
, (b)
$Pe=0.5$
,
$Pe=2.0$
and (d)
$Pe=5.0$
. The dashed lines denote the solid–liquid interface. Qualitatively similar behaviour is also seen for larger values of
$d_{0}$
.

The results discussed here should be contrasted with those for mushy layers, where stability of the system is a non-monotonic function of the strength of the external shear flow, because of the induced flow within the mushy layer (Neufeld & Wettlaufer Reference Neufeld and Wettlaufer2008a
,Reference Neufeld and Wettlaufer
b
). Here, there is no such induced flow and the shear flow only damps perturbations. Namely, the destabilizing factor here is the convective flow that tends to melt the solid phase by enhanced heat transport. The effect of the shear flow is to reduce the strength of vertical motions and hence the upward heat transport. This leads to the decay of perturbations for
$Pe\geqslant 0.22$
.
5 Conclusions
We have studied the effects of shear and buoyancy-driven flow of a pure melt over its solid phase. A linear stability analysis shows that buoyancy is the only destabilizing factor in the system. Shear flow stabilizes the system by reducing the strength of vertical motions and hence vertical heat transport by the convective flow. Our calculations show that for
$Pe$
as small as
$0.22$
, all modes of perturbation decay and the growth rate asymptotes to a negative value for
$Pe\geqslant 1$
. However, we point out the interesting possibility of a re-entrant interfacial instability at much larger
$Pe$
, where shear flow instabilities coupled with buoyancy might enhance vertical motions.
There are clearly implications for situations in which there is a shear flow over a dissolving phase boundary accompanied by a temperature gradient, so that there are potentially three interacting fields (momentum, compositional and thermal) of influence. Pressure fluctuations associated with interfacial corrugations in mushy layers exposed to shear flow can be relieved by dissolution and solidification of the mushy layer itself (Neufeld & Wettlaufer Reference Neufeld and Wettlaufer2008a ,Reference Neufeld and Wettlaufer b ). However, when the solid phase is pure, as in the case studied here, imposing a shear flow with impurities and superheat should lead to interesting phenomena since the temperature of maximum density of aqueous solutions depends on impurity concentration.
Acknowledgements
The authors acknowledge the support of the University of Oxford and Yale University. S.T. acknowledges a Research Fellowship from All Souls College, Oxford, a NASA Graduate Research Fellowship and helpful discussions with A. J. Wells. J.S.W. acknowledges NASA Grant NNH13ZDA001N-CRYO, Swedish Research Council grant no. 638-2013-9243, and a Royal Society Wolfson Research Merit Award for support.
Appendix A. Re-analysis of the experimental results of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980)
In this section, we reanalyse the experimental data of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) to show the effects of buoyancy on their observations of interfacial instability. We adopt the same notation.
A.1 Details of the experiments
Figure 11 shows a schematic of the experimental study of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980). The bottom wall is maintained at a temperature
$T_{w}$
, the ice–water interface is at the bulk equilibrium temperature
$T_{f}$
, and the far-field temperature,
$T_{\infty }$
, is such that
$T_{\infty }>T_{f}>T_{w}$
. The far-field flow speed is
$U_{\infty }$
. Because pure water has a density maximum at
$4\,^{\circ }$
C, the water column above the ice layer is unstably stratified.
A.2 Monin–Obukhov (M–O) theory for a smooth surface
In wall-bounded turbulent shear flows of neutrally buoyant fluids, the flow consists of the inner and the outer regions (Monin & Yaglom Reference Monin and Yaglom1971; Sreenivasan Reference Sreenivasan1989). The inner region is subdivided into: (i) the viscous sublayer, which is closest to the wall, where the effects of viscosity are dominant; (ii) the buffer layer, which is next to the viscous sublayer, where viscous and inertial effects are equally important; and (iii) the log layer, where neither the effects of the wall nor that of the outer region are important. In the limit of asymptotically large Reynolds number, scaling arguments for the behaviour of the mean horizontal velocity,
$U$
, in the log layer lead to (Monin & Yaglom Reference Monin and Yaglom1971; Sreenivasan Reference Sreenivasan1989)
where
$k_{s}=0.41$
is the Kármán constant, and
$B=5.5$
is another constant. The constants
$k_{s}$
and
$B$
are believed to be universal, but their values have been determined only empirically (Sreenivasan Reference Sreenivasan1989). The superscript
$^{+}$
denotes non-dimensionalization by
$u_{\ast }$
, the friction velocity, and
$l_{v}=\unicode[STIX]{x1D708}/u_{\ast }$
, the viscous length scale.
Schematic of the experimental study of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980).

In the case of wall-bounded shear flows of stratified fluids, the stratification affects the mean velocity as follows. If the flow is unstably stratified, there are more vigorous vertical motions and thus more vertical mixing. Hence, the mean velocity at any location is smaller than that for a neutrally buoyant fluid at the same location. However, if the flow is stably stratified then vertical motions are suppressed, leading to a mean velocity that is larger than that for a neutrally buoyant fluid (e.g. Monin & Yaglom Reference Monin and Yaglom1971; Turner Reference Turner1979).
Velocity measurements from figure 4 of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) for
$Re_{\unicode[STIX]{x1D6FF}}=11\,000$
, where
$Re_{\unicode[STIX]{x1D6FF}}$
is the Reynolds number based on the thickness of the boundary layer. From left to right, the positions in the inset are numbered 1, 2, 3 and 4. The solid line is the law of the wall (Monin & Yaglom Reference Monin and Yaglom1971).

Figure 12 shows
$U^{+}(y^{+})$
at different locations in the experiments of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980); the mean velocity profiles show a systematic deviation from the log layer, indicating unstable stratification. Moreover, the amplitude of the rippled interface they observed was small compared to its wavelength. Hence, we treat the surface as planar for the purpose of quantifying the effects of stratification, for which we extend the M–O theory.
The relative effects of inertia and buoyancy are represented by the M–O length scale, denoted by
${\mathcal{L}}$
(Monin & Yaglom Reference Monin and Yaglom1971). For stable stratification
${\mathcal{L}}>0$
and for unstable stratification
${\mathcal{L}}<0$
, with the effects of stratification being important for distances
$y>O(|{\mathcal{L}}|)$
from the wall. Following Monin & Yaglom (Reference Monin and Yaglom1971), we let
$\unicode[STIX]{x1D709}=y/{\mathcal{L}}$
and write
where
$f$
is an unknown function of
$\unicode[STIX]{x1D709}$
and
$\unicode[STIX]{x1D719}(\unicode[STIX]{x1D709})=\unicode[STIX]{x1D709}f(\unicode[STIX]{x1D709})$
. Scaling equation (A 2) with
$u_{\ast }$
and
$l_{v}$
, we have
For
$y^{+}/{\mathcal{L}}^{+}\ll 1$
,
$\unicode[STIX]{x1D719}$
can be expanded in a power series:
$\unicode[STIX]{x1D719}=1+\unicode[STIX]{x1D6FD}(y^{+}/{\mathcal{L}}^{+})+$
higher order terms. Using this in equation (A 3) and integrating with respect to
$y^{+}$
gives
In the limit
${\mathcal{L}}\rightarrow \infty$
, equation (A 4) should reduce to the classical law of the wall, which gives
$A=B=5.5$
. As this analysis is valid for distances ‘far away’ from the wall, the value of
$\unicode[STIX]{x1D6FD}$
is taken to be
$0.6$
(Monin & Yaglom Reference Monin and Yaglom1971).
A.3 Comparison with the experiments
The velocity profile given by equation (A 4) can now be fit to the data of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980). Here,
$b=\unicode[STIX]{x1D6FD}/(k_{s}{\mathcal{L}}^{+})$
is the only fitting parameter, with
$\unicode[STIX]{x1D6FD}$
and
$k_{s}$
already known. Figure 13 shows the fits of (A 4) to the data in figure 12.
Comparison of the theory (A 4) with the measurements of Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) at (a) position 1, (b) position 2, (c) position 3 and (d) position 4. Circles: data from Gilpin et al. (Reference Gilpin, Hirata and Cheng1980); dashed line:
$U^{+}=1/k_{s}\log (y^{+})+A$
; solid line:
$U^{+}=1/k_{s}\log (y^{+})+by^{+}+A$
. Here,
$b=-0.00092$
,
$-0.0015$
,
$-0.0012$
and
$-0.0008$
at positions 1, 2, 3 and 4, respectively.

The averaged value of
$b$
from the fits to the data at the four positions is
$b_{avg}=-0.0011$
, and hence
${\mathcal{L}}_{avg}^{+}=\unicode[STIX]{x1D6FD}/(k_{s}b)=-1330.38$
. Thus, because
${\mathcal{L}}_{avg}<0$
, we confirm that the water column was unstably stratified. From the range of values given for the free-stream velocity in the Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) experiments, we take
$U_{\infty }=0.5~\text{ms}^{-1}$
and use their equation (13),
to obtain
$u_{\ast }=0.033~\text{ms}^{-1}$
for
$Re_{\unicode[STIX]{x1D6FF}}=11\,000$
. Taking
$\unicode[STIX]{x1D708}=10^{-6}~\text{m}^{2}\text{s}^{-1}$
, we obtain
$l_{v}=29.93~\unicode[STIX]{x03BC}\text{m}$
, and hence
$|{\mathcal{L}}_{avg}|=|{\mathcal{L}}_{avg}^{+}|l_{v}=0.04$
m. The height of the test section reported is
$0.457$
m, which makes
$|{\mathcal{L}}_{avg}|$
about
$9\,\%$
of the test-section height. However, because of the departure of velocity profiles from the classical log law at smaller distances from the ice surface, the value of
$|{\mathcal{L}}_{avg}|$
estimated here may be larger than the actual value.
Finally, because the
$Re$
and (estimated)
$Ra$
for their experiments indicate that their flow was in a turbulent regime, our calculations lead us to speculate that the instability observed by Gilpin et al. (Reference Gilpin, Hirata and Cheng1980) may be ostensibly nonlinear; a topic for future study.













