1. Introduction
Melting and freezing processes between ice and water play an important role in the environment. For instance, the melting of ice shelves, i.e. floating glacial ice, can lead to reduced buttressing of the grounded polar ice sheets and increased sea-level rise (Pritchard et al. Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van Den Broeke and Padman2012; Rignot et al. Reference Rignot, Jacobs, Mouginot and Scheuchl2013; Kennicutt Reference Kennicutt2019), while freezing of high-latitude oceans by a cold atmosphere results in sea-ice formation, increased albedo and increased ocean salinity through brine rejection (Wells, Hitchen & Parkinson Reference Wells, Hitchen and Parkinson2019). Icebergs, ice shelves and sea ice are kilometre-scale objects with long lifetimes but their evolution is controlled by heat and salt fluxes across millimetre-thin ice–water boundary layers, which fluctuate rapidly (Dinniman et al. Reference Dinniman, Asay-Davis, Galton-Fenzi, Holland, Jenkins and Timmermann2016). The front of a marine-terminating glacier can melt as fast as several metres per day horizontally (as recently reported for the LeConte Glacier, Sutherland et al. Reference Sutherland, Jackson, Kienholz, Amundson, Dryer, Duncan, Eidam, Motyka and Nash2019), but an ice shelf around Antarctica typically melts at a rate of only a few centimetres per day or less (Dutrieux et al. Reference Dutrieux, Stewart, Jenkins, Nicholls, Corr, Rignot and Steffen2014). On the other hand, ocean currents are most often turbulent and exhibit temporal variabilities on fast time scales of just a few seconds (Davis & Nicholls Reference Davis and Nicholls2019), such that phase changes between ice and water are multi-physics phenomena with large scale separation.
An important consequence of phase changes is that topographical features can emerge at the ice–water interface when the rate of melting and freezing is spatially variable. Basal channels (Stanton et al. Reference Stanton, Shaw, Truffer, Corr, Peters, Riverman, Bindschadler, Holland and Anandakrishnan2013; Gourmelen et al. Reference Gourmelen2017) and terraces (Dutrieux et al. Reference Dutrieux, Stewart, Jenkins, Nicholls, Corr, Rignot and Steffen2014) have been observed at hundreds-of-metre to kilometre scales under ice shelves, the underside of icebergs exhibit ablation channels at the metre scale and scallops at the tens-of-centimetre scale (Hobson, Sherman & McGill Reference Hobson, Sherman and McGill2011), and, more generally, rough features can be seen from the centimetre scale to tens-of-metre scale under sea ice (Wadhams, Wilkinson & McPhail Reference Wadhams, Wilkinson and McPhail2006; McPhee Reference McPhee2008; Lucieer et al. Reference Lucieer, Nau, Forrest and Hawes2016) and up to the kilometre scale under ice shelves (Nicholls Reference Nicholls2006). The interplay between flow dynamics and phase changes leading to the generation and persistence of topographical features in the environment is of fundamental importance. The presence of topography can significantly affect the long-term flow dynamics as well as the average melting or freezing rate of the ice boundary, as suggested by, for example, the large spatial variability of melting of basal terraces (Dutrieux et al. Reference Dutrieux, Stewart, Jenkins, Nicholls, Corr, Rignot and Steffen2014) and recent laboratory experiments on ice scallops (Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019).
Buoyancy forces play an important role in the coupling between phase changes, flow dynamics and topography generation. Buoyancy forces can be stabilizing or destabilizing depending on the relative orientation between the gravitational acceleration vector and the density gradient. In polar seas, cold and fresh melt water near the ice boundary is lighter than the surrounding water, such that buoyancy forces are restoring below a horizontal ice boundary (e.g. below an ice shelf) and drive upwellings along a vertical ice face (e.g. at the front of a marine-terminating glacier). In a cold freshwater system the thermal expansion coefficient of water is negative, i.e. it is negative for temperatures 
$0\,^\circ \textrm {C} < T < 4\,^\circ$C at atmospheric pressure (Thoma et al. Reference Thoma, Grosfeld, Smith and Mayer2010), such that buoyancy forces are stabilizing for water under an ice cover (e.g. as in a frozen lake) but destabilizing for water above ice (e.g. for a supraglacial lake or river).
For a system dominated by destabilizing buoyancy forces, the interplay is strong between fluid dynamics and phase topography. In a thermally stratified fluid with finite depth below a solid phase, the unstable density stratification sets up a large-scale circulation known as Rayleigh–Bénard convection with alternating warm upwelling and cold downwelling regions. The warm upwellings drive stronger melting than the cold downwellings, such that a topography can emerge from spatially variable heat fluxes. The topography enhances the large-scale circulation, such that a positive feedback is obtained and the flow dynamics and solid boundary become phase locked (Rabbanipour Esfahani et al. Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018; Favier, Purseed & Duchemin Reference Favier, Purseed and Duchemin2019). Dissolution of a phase boundary, i.e. with phase changes driven by concentration gradients rather than temperature effects, in a gravitationally unstable fluid can also lead to convective motions and the generation of three-dimensional topography (even in the absence of a large-scale circulation), as shown by experiments (Kerr Reference Kerr1994; Sullivan, Liu & Ecke Reference Sullivan, Liu and Ecke1996) and numerical simulations (Philippi et al. Reference Philippi, Berhanu, Derr and du Pont2019). Streamwise patterns also emerge in dissolution experiments when the phase boundary is not perpendicular to gravity but inclined (Allen Reference Allen1971; Cohen et al. Reference Cohen, Berhanu, Derr and du Pont2020).
Despite the existence of many studies on pressure-driven and shear flows (Kelly Reference Kelly1994; Zonta & Soldati Reference Zonta and Soldati2018), the possibility for topography to emerge between a horizontal boundary layer flow and a solid phase, i.e. perpendicular to gravity, is not well understood, at least compared with the case of topography generation by Rayleigh–Bénard convection. Boundary layer flows strongly affected by shear, such as under ice shelves, are yet as common (if not more) as buoyancy-driven flows in the environment, such that there is significant interest in predicting their ability to generate topographical features (or roughness) at horizontal ice boundaries and the impact the sustained topography can have on overall melt rates. Using laboratory experiments, Gilpin, Hirata & Cheng (Reference Gilpin, Hirata and Cheng1980) demonstrated the existence of an interfacial instability and the generation of ripples at an ice boundary below a horizontal turbulent boundary layer water flow. The experiments had a modest unstable density stratification owing to the negative thermal expansion coefficient of freshwater at low temperatures (Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019). However, Bushuk et al. (Reference Bushuk, Holland, Stanton, Stern and Gray2019), who conducted an experiment similar to Gilpin et al. (Reference Gilpin, Hirata and Cheng1980), argued that the buoyancy forcing can be neglected in the large velocity regime of the experiments, thus confirming the possibility for topography generation in the absence of vertical convection. The necessary condition for an interfacial instability to develop, regardless of the type of density stratification, is that the maximum of mass transfer from the solid to the fluid due to a heat flux or concentration gradient (resulting in ablation) at the boundary be shifted by 
$-90^\circ$ to 
$+90^\circ$ compared with the maximum (crest) of boundary topography (Hanratty Reference Hanratty1981). Recently, Claudin, Durán & Andreotti (Reference Claudin, Durán and Andreotti2017) demonstrated that such a shift is possible for a horizontal neutral turbulent flow and proposed a saturation mechanism for the finite amplitude of two-dimensional scallops. Three-dimensional effects and buoyancy forces are expected to play an important role in topography generation and melting rates but were not considered in the study of Claudin et al. (Reference Claudin, Durán and Andreotti2017), which also relied on parameterized flow nonlinearities. Thus, additional efforts are necessary to improve our understanding of the physical mechanisms leading to the generation and saturation of phase topography by horizontal shear flows.
Here, we investigate the generation of topography at a phase boundary adjacent to a shear flow affected by buoyancy via direct numerical simulations. We focus on the case of an initially flat and horizontal solid, i.e. perpendicular to gravity, and investigate the influence of density stratification on the topography obtained and the coupled fluid–solid dynamics. Our numerical model solves the evolution of the fluid and solid phases simultaneously using the phase-field method. The phase-field method is a one-domain two-phase fixed-grid method that was originally developed by the metallurgy community for relatively smooth flows (Wang et al. Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993; Karma & Rappel Reference Karma and Rappel1998; Beckermann et al. Reference Beckermann, Diepers, Steinbach, Karma and Tong1999), but which was recently applied to the case of vigorous convective flows (Favier et al. Reference Favier, Purseed and Duchemin2019; Purseed et al. Reference Purseed, Favier, Duchemin and Hester2020). Other methods that simultaneously solve for the evolution of a fluid phase and a solid phase include the enthalpy method (Ulvrová et al. Reference Ulvrová, Labrosse, Coltice, Råback and Tackley2012), the level-set method (Gibou et al. Reference Gibou, Chen, Nguyen and Banerjee2007), the lattice-Boltzmann method (Rabbanipour Esfahani et al. Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018) and two-domain moving-boundary methods (Ulvrová et al. Reference Ulvrová, Labrosse, Coltice, Råback and Tackley2012). The main advantage of the phase-field method over these other methods is that it can be implemented relatively easily in any fluid solver.
Our study aims to contribute to the physical understanding of topography generation by shear flows at horizontal boundaries and the associated changes in mean melt rates, as most recently investigated theoretically by Claudin et al. (Reference Claudin, Durán and Andreotti2017) and experimentally by Bushuk et al. (Reference Bushuk, Holland, Stanton, Stern and Gray2019). Numerical constraints force us to consider an idealized set-up however, such that our fluid and solid phases are not exactly representative of water and ice. Notably, we assume that the fluid and solid have the same thermodynamical and transport properties, e.g. the same thermal conductivity, and we consider an anomalously warm fluid in order to minimize the time scale separation between the turbulent dynamics and generation of boundary topography. Due to computational constraints, the external flow in our simulations is also weaker than what may be expected for scallops to emerge (Claudin et al. Reference Claudin, Durán and Andreotti2017; Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019).
The main result of our paper is that topographical features spontaneously emerge at the ice–water interface due to uneven melting of the solid boundary by the shear flow. We investigate the effect of background density stratification and we demonstrate that the topography is dominated by keels and channels that are aligned with the direction of the mean flow in all cases.
We organize the manuscript as follows. In § 2 we describe the phase-field method, the dimensionless equations and the numerical method. In § 3 we present and discuss the direct numerical simulation results obtained for three different background stratifications. In § 4 we discuss the link between our results and geophysical applications and explain why we did not observe three-dimensional scallops. In § 5 we conclude. Finally, in appendices A–D, we provide additional details about the method and results.
2. Model
2.1. Phase-field method
We investigate the generation of topography due to uneven melting and freezing at a fluid–solid interface. The solid is fixed and located above the fluid where a Poiseuille/channel flow develops due to an external pressure gradient (see figure 1). The initial thickness of the fluid (respectively solid) layer is 
$H$ (respectively 
$H/2$), such that the channel full depth is 
$3H/2$. The domain length (in the direction of the flow) is 
$L_x=4{\rm \pi} H$ and the transverse width is 
$L_y=2{\rm \pi} H$. We define a Cartesian coordinate system 
$(x,y,z)$ centred on the bottom of the channel with the 
$z$-axis vertically upward, i.e. opposite to gravity, and use superscripts 
$^{(\,f)}$ and 
$^{(s)}$ to denote variables in the fluid and the solid, respectively. The fluid velocity 
${\boldsymbol {u}}^{(\,f)}$ and pressure 
$p^{(\,f)}$ evolve according to the Navier–Stokes equations under the Boussinesq approximation. For simplicity, we assume that the solid and fluid phases have the same thermodynamical and transport properties, i.e. the same reference density 
$\rho _f$, the same specific heat capacity at constant pressure 
$c_p$ and the same thermal conductivity 
$k$. Thus, the temperatures 
$T^{(\,f)}$ and 
$T^{(s)}$ evolve according to the same advection–diffusion (energy) equation, which turns into the heat equation where the velocity is zero. We consider a generic linear equation of state for the fluid, i.e. not specific to water, with the density related to temperature through
\begin{equation} \rho^{(\,f)}=\rho_f\left(1-\alpha T^{(\,f)}\right), \end{equation}
with 
$\alpha$ the thermal expansion coefficient. For a pure component flow, the fluid–solid interface must be at the temperature of melting, denoted by 
$T_m$, and the movement of the interface is governed by the Stefan condition, i.e.
\begin{gather} T^{(\,f)}=T^{(s)}=T_m, \end{gather}
\begin{gather}\rho_s\mathcal{L}v_n= q_n^{(\,f)}-q_n^{(s)} = -k_f\hat{\textbf{n}}\cdot\boldsymbol{\nabla} T^{(\,f)}+k_s\hat{\textbf{n}}\cdot\boldsymbol{\nabla} T^{(s)}, \end{gather}
where 
$\rho _s$ is the reference density of the solid, 
$v_n$ is the interface velocity in the direction normal to the interface and directed toward the solid phase (supported by unit vector 
$\hat {\boldsymbol {n}}$), 
$\mathcal {L}$ is the latent heat of fusion per unit mass, 
$q_n$ is the heat flux in direction 
$\hat {\boldsymbol {n}}$, 
$k_s$ (respectively 
$k_f$) is the thermal conductivity of the solid (respectively fluid) and 
$\boldsymbol {\nabla }$ is the gradient operator (Worster Reference Worster2000). We recall that we assume the same properties for the two phases, i.e. in our case 
$k_f=k_s=k$ and 
$\rho _s=\rho _f$ in (2.2). Note that the properties of water and ice are different under typical atmospheric pressure and near-freezing temperature conditions, i.e. 
$\rho _f\approx 999$ kg m
$^{-3}$, 
$c_{pf}\approx 4200$ J kg
$^{-1}$ K
$^{-1}$ and 
$k_f\approx 0.6$ W m
$^{-1}$ K
$^{-1}$, while 
$\rho _s\approx 917$ kg m
$^{-3}$, 
$c_{ps}\approx 2100$ J kg
$^{-1}$ K
$^{-1}$ and 
$k_s\approx 2.2$ W m
$^{-1}$ K
$^{-1}$. The relative differences are small however, i.e. within a factor of four or less, such that we do not expect fundamental differences between the physics in our model and real processes involving water and ice in nature.
(a) Simulation snapshot showing the temperature field in the fluid (red colourmap) and the solid (blue colourmap) for neutral stratification (
$Ri_*=0$) at late time, i.e. such that it is representative of the statistical steady state. The arrows along the four vertical transects and across the front face 
$x=l_x$ (displayed over a limited number of points for clarity) show the velocity vectors. The velocity vectors all point primarily in the direction of the mean flow, 
$+\hat {x}$. (b) Zoom-in on a region of (a). (c) Variation of the phase field from 
$\phi =1$ in the fluid to 
$\phi =0$ in the solid along the vertical solid line drawn in (b). The non-dimensional lengths in the 
$x$, 
$y$ and 
$z$ directions are 
$l_x=4{\rm \pi}$, 
$l_y=2{\rm \pi}$ and 
$l_z=1.5$, respectively.
Here we use a volume-penalization method (Angot, Bruneau & Fabrie Reference Angot, Bruneau and Fabrie1999), which is a type of immersed boundary method, combined with the phase-field method, in order to solve for phase-change processes and the evolution of the variables in the fluid and the solid simultaneously. Specifically, we solve the Navier–Stokes equations in the Boussinesq approximation and the advection–diffusion (energy) equation for temperature combined with an equation for the fluid fraction 
$\phi$, i.e.
\begin{gather} \frac{\partial {\boldsymbol{u}}}{\partial t} + \phi \left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} \right){\boldsymbol{u}} = \nu \nabla^2{\boldsymbol{u}} - \frac{1}{\rho_f}\boldsymbol{\nabla} p + \alpha g T \hat{\textbf{z}} + \frac{\varPi}{\rho_f}\hat{\textbf{x}} - \frac{(1-\phi)}{\tau_p}{\boldsymbol{u}}, \end{gather}
\begin{gather}\frac{\partial T}{\partial t} + \phi \left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} \right)T = \kappa \nabla^2T - \frac{\mathcal{L}}{c_p}\frac{\partial \phi}{\partial t}, \end{gather}
\begin{gather}\frac{\partial \phi}{\partial t} = a \nabla^2 \phi + b\phi(1-\phi)\left[ 2\phi-1+c(T-T_m)\right], \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}}=0, \end{gather}
with 
${\boldsymbol {u}}=(u,v,w)$, 
$p$ and 
$T$ defined in both the fluid and the solid, i.e. assuming that the fluid and solid phases form a single domain, such that we can drop the superscripts 
$^{(\,f)}$ and 
$^{(s)}$. In (2.3), 
$\nu$ is the constant kinematic viscosity, 
$g$ is the standard gravity, 
$\varPi$ is the imposed pressure-gradient force and 
$\kappa =k/\rho _f/c_p$ is the constant thermal diffusivity; 
$\tau _p$, 
$a$, 
$b$ and 
$c$ are parameters related to volume penalization and the phase-field method, which we define later, and 
$\hat {\boldsymbol {z}}$ and 
$\hat {\boldsymbol {x}}$ are the unit vectors of the 
$z$- and 
$x$-axes, respectively. Note that the third term on the right-hand side of (2.3a) represents the buoyancy force.
The fluid fraction, 
$\phi$, also known as the phase-field variable or order parameter, satisfies a forced diffusion equation (2.3c) with parameters tuned such that 
$\phi$ transitions continuously from 1 in the fluid to 0 in the solid, across a diffuse interface whose thickness is artificial and must be smaller than all physical length scales in the problem, including the viscous length scale (cf. appendix A). The fluid fraction 
$\phi$ is introduced in the momentum (2.3a), energy (2.3b) and continuity (2.3c) equations, in order to modulate locally the importance of each physical process based on the component's phase. For instance, the last term on the right-hand side of (2.3a) is a linear (penalization) damping term, which is active in the solid but inactive in the fluid, while the second term on the right-hand side of the energy equation (2.3b) is a heat sink or source that represents the consumption or release of latent heat associated with melting or freezing. In the limit of an infinitesimally small diffuse interface thickness of the phase field, it has been shown that the dynamics of the fluid–solid interface governed by (2.3) converges to the exact Stefan conditions (2.2), and that the fluid velocity converges to 0 at the fluid–solid interface, thus mimicking a no-slip boundary. Here we multiply by 
$\phi$ the advective terms in (2.3a) and (2.3b), such that they are zero in the solid phase. Previous studies have used both damped and undamped advective terms and we discuss the impact of our choice on the results in appendix B.
2.2. Dimensionless equations
Equations (2.3) can be non-dimensionalized in order to identify the set of independent control parameters. Following previous studies (e.g. Favier et al. Reference Favier, Purseed and Duchemin2019), we use the thermal diffusive time scale 
$\tau _{\kappa }=H^2/\kappa$ as a normalizing time scale, i.e. the dimensionless variables, denoted by tildes, are defined as
\begin{equation} \left.\begin{aligned} & (x,y,z)=(H\tilde{x},H\tilde{y},H\tilde{z}),\quad t=\tau_{\kappa}\tilde{t},\quad u=u_{\kappa}\tilde{u},\\ & T=T_m+{\rm \Delta} T\tilde{T},\quad p = \rho_f u_{\kappa}^2\tilde{p},\quad \phi=\tilde{\phi}, \end{aligned}\right\} \end{equation}
with 
$u_{\kappa }=\kappa /H$ the diffusion velocity scale, and 
${\rm \Delta} T=T_b-T_m$ is the temperature scale with 
$T_b$ the dimensional temperature on the bottom boundary. The time scale 
$\tau _{\kappa }$ is particularly relevant for discussing the long-term dynamics of the system since temperature evolves in the solid through diffusion. We will use the shorter friction time scale to describe relatively rapid processes, such as convection in the fluid (see § 2.3).
Substituting variables (2.4a–f) into (2.3) and dropping tildes, we obtain the dimensionless equations
\begin{gather} \frac{\partial {\boldsymbol{u}}}{\partial t} + \phi \left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} \right){\boldsymbol{u}} = Pr \nabla^2{\boldsymbol{u}} - \boldsymbol{\nabla} p + Pr Ra T \hat{\textbf{z}} + 2Pr^2 Re \hat{\textbf{x}} - Pr\frac{(1-\phi)}{\varGamma}{\boldsymbol{u}}, \end{gather}
\begin{gather}\frac{\partial T}{\partial t} + \phi \left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} \right)T = \nabla^2T - St\frac{\partial \phi}{\partial t}, \end{gather}
\begin{gather}\frac{\partial \phi}{\partial t} = A \nabla^2 \phi + B\phi(1-\phi)(2\phi-1+CT), \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}}=0. \end{gather}
The control parameters in (2.5) are the Prandtl number, 
$Pr$, which compares kinematic viscosity to thermal diffusivity, the centreline Reynolds number, 
$Re$, which compares the pressure-gradient force to viscous dissipation, the Rayleigh number, 
$Ra$, which compares buoyancy forces to viscous and thermal dissipation, and the Stefan number, 
$St$, which compares the available sensible heat to the latent heat. They are related to the physical parameters through
\begin{equation} Pr=\frac{\nu}{\kappa}, \quad Re=\frac{\varPi H^3}{2\rho_f\nu^2},\quad Ra = \frac{\alpha g {\rm \Delta} T H^3}{\nu\kappa},\quad St=\frac{\mathcal{L}}{c_p{\rm \Delta} T}. \end{equation}
The additional parameters 
$\varGamma =\tau _p\nu H^2/\kappa ^2$, 
$A=a/\kappa$, 
$B=b/(\kappa /H^2)$ and 
$C$ are non-physical and prescribed based on numerical constraints of the volume-penalization and phase-field methods (cf. appendix A). The problem is fully specified once 
$Pr$, 
$Re$, 
$Ra$ and 
$St$ are known and the boundary conditions are prescribed. Here we enforce a no-slip, fixed-temperature condition at the top of the ice, i.e. 
${\boldsymbol {u}}=\boldsymbol {0}$ and 
$T=T_t<0$ at 
$z=1.5$. We impose free-slip, fixed-temperature conditions on the bottom boundary, i.e. 
$\partial _zu=\partial _zv=w=0$ and 
$T=1$ at 
$z=0$, such that we simulate only one half of a full channel flow (to reduce computational costs). The dimensionless melting temperature is 
$T=0$. The initial interface position is 
$z=1$ and we note that 
$(l_x,l_y,l_z)=(4{\rm \pi} ,2{\rm \pi} ,1.5)$, the domain lengths in dimensionless space. The initial condition in the fluid is a half-channel laminar Poiseuille flow superimposed with divergence-free white noise for the velocity fluctuations.
We will generally discuss our results in terms of the steady-state friction (or shear) Reynolds number, 
$Re_*$, and the friction Richardson number, 
$Ri_*$, i.e.
\begin{equation} Re_*=\sqrt{2Re}=\sqrt{\frac{\varPi H^3}{\rho_f\nu^2}},\quad Ri_*=\frac{-Ra}{PrRe}=\frac{-2\rho_f\alpha{\rm \Delta} T g}{\varPi}, \end{equation}
since they are more commonly used than 
$Re$ and 
$Ra$ in turbulent channel flow studies (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Zonta & Soldati Reference Zonta and Soldati2018). The key difference between 
$Re$ and 
$Re_*$ is that the former is based on the velocity on the bottom free-slip boundary of the channel in the laminar regime, which is 
$[\varPi H^2(1-z^2/H^2)/(2\rho _f\nu )]|_{z=0}$ using dimensional variables, while the latter is based on the friction velocity, which is 
$\sqrt {-\tau _w}=\sqrt {\varPi H/\rho _f}$, with 
$\tau _w$ the mean wall shear stress, again using dimensional variables. Here, we favour the friction Richardson number over the Rayleigh number as the control parameter, even when the stratification is unstable, because they are both input parameters and because the wall shear stress is an important driver of turbulence in all cases. The importance of shear forces compared with buoyancy forces can be estimated from the Monin–Obukhov length, which is
\begin{equation} L_{MO} = \frac{Re_*^3Pr^2}{Nu\,Ra}, \end{equation}
in terms of dimensionless variables and which is often reported in mixed-convection experiments (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020), with 
$Nu$ the Nusselt number, which we define later (see (2.9a–c)). The Monin–Obukhov length estimates the distance from the boundary within which shear is as important or more important than buoyancy. In our simulations, we always have 
$L_{MO}>0.97$, such that shear is a significant source of turbulence throughout the domain.
We investigate the effect of background density stratification on the generation of topography at the fluid–solid interface by considering three distinct values of 
$Ri_*$, i.e. 
$Ri_*=40$, 
$Ri_*=0$ and 
$Ri_*=-40$, for which the stratification is stable, neutral and unstable, respectively. For simplicity and computational expediency, all other parameters (except 
$T_t$) are fixed such that the flow is (moderately) turbulent and phase changes are relatively rapid, i.e. we set 
$Re_*=150$ (
$Re=11\,250$), 
$Pr=1$ and 
$St=1$. For each 
$Ri_*$, we set 
$T_t$ such that the initial heat flux in the ice, 
$-T_t/2$, is almost equal to the heat flux in the fluid when there is no melting. As a result, the fluid–solid interface position does not move significantly in time (at least initially) and we can maximize numerical resolution around the interface with a fixed grid. For reference, the Rayleigh number for the unstable stratification case (
$Ri_*=-40$) is 
$Ra=4.5\times 10^5$, which is above the instability onset for Rayleigh–Bénard convection rolls in the streamwise direction (
$Ra\approx 1101$) of thermally stratified plane Poiseuille flow (Chandrasekhar Reference Chandrasekhar1961; Gage & Reid Reference Gage and Reid1968). Note that changing the sign of 
$Ri_*$ can be obtained by changing the sign of the thermal expansion coefficient 
$\alpha$, which indeed can be either positive or negative depending on the fluid state. We discuss the geophysical relevance of our choice of parameters, including 
$\alpha$, in § 4.
We solve (2.5) using the open-source pseudo-spectral direct numerical simulation (DNS) code DEDALUS (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). We use 256 Fourier modes in the 
$x$ and 
$y$ directions and a compound Chebyshev basis with a total of 288 modes in the 
$z$ direction unless stated otherwise. The use of a compound Chebyshev basis allows us to have a stretched grid in the vertical direction with refined (respectively coarse) resolution near the mean fluid–solid interface (respectively in the fluid bulk). Here, the Chebyshev collocation grid has a resolution equal to approximately one-fourth of a wall unit 
${\rm \Delta} z^+=1/Re_*\approx 0.0066$ at and around the fluid–solid interface and equal to approximately one wall unit in the fluid bulk, whereas the Fourier collocation grids have a uniform resolution of roughly 7 and 3.5 wall units in 
$x$ and 
$y$, which is within the recommended resolution for channel flow simulations (see e.g. Moin & Mahesh (Reference Moin and Mahesh1998) and appendices A and C for more details). We use a two-step implicit–explicit Runge–Kutta scheme for time integration. The Courant–Friedrichs–Lewy condition is typically set to 0.2 in the transient initial stage and 0.4 later on. At statistical steady state, the time step is typically 
$10^3\text {--}10^4$ times smaller than the friction time scale 
$1/(Re_*Pr)$, which is equal (in terms of dimensional variables) to 
$H$ divided by the steady-state friction velocity. We run each simulation for approximately 4 diffusive time scales, or 600 friction time scales, which takes roughly 2 million time steps, such that the total cost of the simulations is of the order of 1 million CPU hours. Figure 1(a) shows a snapshot of the temperature field in the fluid (red colourmap) and the solid (blue colourmap), as well as the velocity vectors (arrows) at select locations for 
$Ri_*=0$. Figure 1(c) shows the variations of the phase field along the thick solid line drawn in figure 1(b). The transition from 
$\phi =1$ in the fluid to 
$\phi =0$ in the solid occurs over a very thin diffuse interface of thickness 
${\approx }0.007$. Simulation parameters and output variables are provided in table 1.
Simulation parameters and selected output variables. The output bulk Reynolds numbers 
$Re_b^I$ and 
$Re_b^{III}$ and Nusselt numbers 
$Nu_b^I$ and 
$Nu_b^{III}$ are volume-averaged and time-averaged over 50 friction time units before the end of stages I and III, respectively. Here 
$q^s$ is the conductive heat flux through the ice imposed as an initial condition at the beginning of stage II; 
$\xi ^{III}$, 
$\xi _-$ and 
$\xi _+$ are the mean interface position, the maximum amplitude of the keels and the maximum amplitude of the channels averaged over 50 friction time units at the end of stage III; 
$C_D^I$ and 
$C_D^{III}$ are the drag coefficients averaged over 50 friction time units at the end of stage I and stage III, respectively. Note that 
$Re_*=150$, 
$Pr=1$ and 
$St=1$ in all simulations. Also, 
$Ri_*=-40$ corresponds to 
$Ra=4.5\times 10^5$.
2.3. Variables of interest
We define the friction velocity, the bulk velocity and the Nusselt number as
\begin{equation} u_* = \langle \phi \rangle\sqrt{-\overline{\left(\tau_d+\tau_{\nu}+\tau_w\right)}}\Biggl|_{z=1.5}, \quad u_b = \frac{\langle u \rangle}{\langle \phi \rangle}, \quad Nu = \langle q \rangle, \end{equation}
respectively (cf. details in appendix B), where the overbar denotes the horizontal average and 
$\langle \cdot \rangle \equiv \int _{\mathcal {V}}\textrm {d}\mathcal {V}/\mathcal {V}$ denotes the volume average, such that 
$\langle \phi \rangle$ is the mean fluid fraction and 
$u_b$ is the bulk velocity of the fluid phase. In (2.9a–c), 
$\tau _d$, 
$\tau _{\nu }$ and 
$\tau _w$ are the linear damping, viscous and Reynolds shear stresses, and 
$q=wT-\partial _zT$ is the heat flux. At statistical steady state, the full shear stress 
$\bar {\tau }=\overline {\left (\tau _d+\tau _{\nu }+\tau _w\right )}$ is approximately a linear function of 
$z$ and 
$q$ is approximately depth invariant, in agreement with channel flow simulations of a pure fluid (cf. appendix B for details on stresses and depth-independent variables using the phase-field method). Since 
$\bar {\tau }$ is a linear function of 
$z$, 
$u_*$ can be estimated from the full shear stress as 
$\sqrt {-\bar {\tau }}$ at any depth as long as it is properly rescaled by the height at which it is estimated. Here, we use the shear stress at the top boundary 
$z=1.5$ in (2.9a–c) for convenience but with pre-multiplying coefficient 
$\langle \phi \rangle \leq 1$, such that 
$u_*$ is truly the friction velocity at the mean interface position (cf. (2.9a–c)). We denote by 
$\xi$ the fluid–solid interface position, where
\begin{equation} \xi(x,y) = \int_0^{l_z} \phi\,\textrm{d} z, \end{equation}
such that 
$\bar {\xi }=\langle \phi \rangle$ (note that one could alternatively define 
$\xi$ such that it satisfies 
$\phi (z=\xi )=0.5$ or 
$T(z=\xi )=0$), and we denote the melt rate by 
$\dot {m}=\partial _t\xi$. The drag coefficient of the fluid–solid boundary is defined as the ratio of the dimensionless wall shear stress 
$u_*^2$ divided by the dynamic pressure 
$u_b^2/2$, i.e.
\begin{equation} C_D=2\left(\frac{u_*}{u_b}\right)^2 \end{equation}
(García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). The temporal fluctuations of the variables of interest will be mainly reported in terms of the friction time 
$t_*=Re_*Prt$. Occasionally, we will show vertical profiles of variables in terms of the distance from the interface, which we denote by 
$\chi (t,x,y)=\xi (t,x,y)-z$.
3. Results
The key findings of our work are that (i) streamwise topographical features emerge from uneven melting and freezing at a phase boundary when the flow is driven by a pressure gradient, and that (ii) the type of density stratification affects the characteristic amplitude and spanwise wavelength of the streamwise patterns. Thus, after a discussion of the evolution of global flow variables in § 3.1, we directly present the results of the topographical features generated at the fluid–solid boundary in § 3.2. We then investigate the interplay between the turbulent flow, the topography and phase changes in §§ 3.3 and 3.4, and finally discuss the evolution of the mean interface position and the statistics of melting and freezing in § 3.5.
3.1. Simulation stages and global flow variables
We show in figure 2 the friction velocity 
$u_*$, the bulk velocity 
$u_b$ and the Nusselt number 
$Nu$ for stable (top figure), neutral (middle figure) and unstable (bottom figure) stratification. Each simulation is broken down into three main stages, which are highlighted by different colours in figure 2 (note that we do not discuss the results shown in grey, which correspond to the spin-up of the fluid phase without buoyancy effects). The first stage of interest (stage I for 
$t_*^{Ib} \leq t \leq t_*^{Ic}$) is shown in blue and corresponds to the spin-up of the fluid phase with buoyancy effects turned on. Importantly, stage I neglects the solid phase, which is substituted with a simple isothermal no-slip boundary, for computational expediency. The second key stage (stage II for 
$t_*^{Ic} < t \leq t_*^{II}$) is shown in orange and corresponds to the part of the simulation that includes the solid phase with volume penalization turned on, but neglects melting or freezing, such that the solid always occupies 
$1\leq z \leq 1.5$ and the phase field is prescribed as 
$\phi =0.5 \left \{ 1 - \tanh \left [ 2(z-1)/\delta \right ] \right \}$, where 
$\delta$ is the thickness of the diffuse interface. The final third stage (stage III for 
$t>t_*^{II}$) is shown in green and highlights results obtained when all effects are considered, i.e. buoyancy is turned on, there is both the fluid and the solid and phase changes are enabled (cf. additional details on the simulation stages in appendix C). The temperature in the solid is initialized at the beginning of stage II as
\begin{equation} T=-q^s (z-1),\quad 1\leq z\leq~1.5, \end{equation}
where 
$q^s$ is the initial conductive heat flux through the solid, by imposing the fixed-temperature condition 
$T=T_t=-q^s/2$ at the top of the solid. The difference between the heat flux in the fluid and the conductive heat flux in the solid in stage II controls whether the fluid–solid interface melts or freezes once phase changes are turned on in stage III. Here, we set 
$q^s$ to be slightly smaller than the heat flux in the fluid at the end of stage I, which we denote by 
$Nu^I$, such that the solid melts slowly at the beginning of stage III in all three simulations (see further discussion in § 3.5). The bulk Reynolds and Nusselt numbers at the end of stages I and III are defined as
\begin{equation} \left.\begin{aligned} & Re_b^I = \int_{t_*^{Ic}-\varDelta_*}^{t_*^{Ic}} \frac{u_b\,\textrm{d} t_*}{Pr \varDelta_*},\quad Re_b^{III} = \int_{t_*^{III}-\varDelta_*}^{t_*^{III}} \frac{u_b\,\textrm{d} t_*}{Pr\varDelta_*},\\ & Nu^I = \int_{t_*^{Ic}-\varDelta_*}^{t_*^{Ic}}\frac{Nu\,\textrm{d} t_*}{\varDelta_*},\quad Nu^{III} = \int_{t_*^{III}-\varDelta_*}^{t_*^{III}} \frac{Nu\,\textrm{d} t_*}{\varDelta_*}, \end{aligned}\right\} \end{equation}
with 
$\varDelta _*=50$, and are reported with 
$q^s$ in table 1. Note that 
$Re_b\ll Re$ because the flow is turbulent and, hence, experiences enhanced friction at the wall compared with the same flow in the laminar regime.
Time history of the friction velocity 
$u_*$ and bulk velocity 
$u_b$ (left axis) and of the Nusselt number 
$Nu$ (right axis) for (a) stable (
$Ri_*=40$), (b) neutral (
$Ri_*=0$) and (c) unstable (
$Ri_*=-40$) density stratification. We report the friction time 
$t_*=Re_*Prt$ on the horizontal axis and show 
$u_b/20$ instead of 
$u_b$ as 
$u_b/20$ and 
$u_*$ have the same order of magnitude. Each full simulation can be broken down into a series of stages of increasing complexity, which are shown by different colours and are separated by vertical dashed lines. The overall procedure is shown schematically in the top panel. There is an intermediate stage 
$t_*\in [t_*^{Ib},t_*^{Ici}]$ in (a) during which 
$Ri_*=20$. See the text and appendix C for more details.
Buoyancy effects are turned off for 
$t_* \leq t_*^{Ib}$, such that the results of figure 2 are exactly the same for all three simulations until 
$t_* = t_*^{Ib}$. Upon turning on buoyancy, i.e. for 
$t_* \geq t_*^{Ib}$ (blue), the Nusselt number and bulk velocity deviate from the neutral case (middle figure), but with opposite behaviours: 
$Nu$ decreases while 
$u_b$ increases with stabilizing buoyancy effects (top figure), and 
$Nu$ increases while 
$u_b$ decreases with destabilizing buoyancy effects (bottom figure). The friction velocity, on the other hand, remains close to 
$u_*/Pr\approx Re_*$ in all three cases. The effect of background stratification on bulk velocity and heat fluxes are well known from channel flow studies (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017), and the important point is that the heat flux is the variable that changes the most with buoyancy effects. Here, 
$Nu^I=3.02$, 4.58 and 6.74 for 
$Ri_*=40$, 0 and 
$-$40, respectively (cf. table 1). It is worth noting that while 
$Nu$ remains the same between stage I and stage II (in a time-average sense), 
$u_*$ and 
$u_b$ show some variations as a result of turning on volume penalization and adding a solid phase. The large dip of 
$u_*$ at 
$t_* \approx t_*^{Ic}$ is merely the result of a sudden deceleration of the mean flow close to the interface, due to the addition of linear damping, which is transient, as can be seen from the rapid return of 
$u_*$ to its statistically steady-state value of 
$u_*\approx 150$. The drop of the bulk velocity is similarly due to the added linear damping. However, unlike the dip in 
$u_*$, the drop in 
$u_b$ persists at all times, implying that volume penalization results in anomalous drag on the mean flow. Here, the relative drop of bulk velocity is of the order of 5 % and the profiles of temperature and velocity close to the fluid–solid interface in stage II reproduce closely those obtained in stage I (see appendix A). Therefore, we consider the discrepancy to be small enough not to warrant a computationally costly increase in resolution or further tuning of the phase-field parameters. When melting is turned on, i.e. for 
$t_*>t_*^{II}$ (green), global flow variables show different behaviours depending on 
$Ri_*$. For the stable case, 
$u_*$, 
$u_b$ and 
$Nu$ exhibit moderately large fluctuations (as in previous stages), but do not exhibit any time-mean deviation (top figure). For the neutral case, we find a small increase in 
$u_*$, 
$u_b$ and 
$Nu$ (middle figure). For the unstable case, we find that 
$u_*$ and 
$u_b$ increase slightly, while 
$Nu$ increases substantially (bottom figure). The analysis presented in the next sections explains these behaviours. Eventually, all simulations reach a statistical steady state.
We show in figure 3 the temporal evolution of another global variable, namely, the drag coefficient, 
$C_D$, which is of significant interest in inferring melt rates from resolved variables in coarse models (using, for instance, the three-equation model; see Holland & Jenkins Reference Holland and Jenkins1999). The drag coefficient is of order 
$10^{-2}$ and decreases (respectively increases) significantly at 
$t_*=t_*^{Ib}$, i.e. when the stratification becomes stable (respectively unstable). The decrease (respectively increase) of 
$C_D$ results from an increase (respectively decrease) of the potential energy barrier in stirring the mean shear and bringing momentum upward with increasing stable (respectively unstable) stratification and is in agreement with previous studies (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). In stage II, 
$C_D$ increases because 
$u_b$ decreases moderately upon turning on volume penalization (cf. figure 2). In stage III, 
$C_D$ has similar values to stages I and II (cf. reported values in table 1), showing that it is not modified by the topographical features obtained in DNS, perhaps because they are aligned with the main flow direction (see § 3.2).
Drag coefficient 
$C_D$ as a function of time 
$t_*$ for unstable, neutral and stable stratification. The colours highlight different simulation stages as in figure 2.
3.2. Spontaneous generation of channels and keels
The mean interface position does not vary significantly in our simulations, due to our choice of initial and boundary conditions for the solid, but uneven melting by the turbulent flow still generates large-amplitude topography, which we discuss in this section. We denote variables averaged in the 
$x$ direction by a tilde (
$\;\widetilde {\;}\;$) and variables averaged in the 
$x$ direction minus the horizontal mean by a prime (
$'$), such that e.g. 
$\xi '=\tilde {\xi }-\bar {\xi }$ represents the spanwise variations of the streamwise-averaged topography around the horizontal mean.
We show snapshots of the two-dimensional fluid–solid interface 
$\xi$ at the end of stage III in figure 4(a–c) for stable, neutral and unstable stratification, respectively (cf. movies 1–3 in the supplementary material available at https://doi.org/10.1017/jfm.2020.1064 to see the temporal evolution of the interface). In all three cases, the topography is dominated by channels (troughs in the solid; brown colour) and keels (excursions of solid into the fluid; green colour) aligned with the streamwise direction. We reach a statistical steady state relatively quickly in all cases after turning on phase changes, such that the patterns in figure 4(a–c) are representative of the interface topography throughout most of stage III (see movies 1–3 in the supplementary material). We show in figure 4(d–f) the Hövmoller diagrams of the channels and keels by plotting 
$\xi '$ in the 
$(t_*,y)$ plane for all of stage III. It can be seen that the characteristic amplitudes of the channels and keels saturate almost immediately for stable and neutral stratification and well before the end of stage III for unstable stratification. The steady-state amplitude of the biggest channels, 
$\xi _+$ (maximum of 
$\xi '$), and the steady-state amplitude of the deepest keels, 
$\xi _-$ (minus the minimum of 
$\xi '$), increase with decreasing 
$Ri_*$ (i.e. from figure (d) to (f)). The crest-to-trough amplitude is roughly 5, 10 and 45 times the viscous length scale 
$\delta _{\nu }=1/Re_*$ for stable, neutral and unstable stratification, respectively (note that 
$\delta _{\nu }$ is roughly equal to the diffuse interface thickness; cf. appendix A). Thus, the crest-to-trough amplitude is of the same order as the viscous sublayer thickness, which is approximately 
$5\delta _{\nu }$, for stable and neutral stratification, but extends beyond the buffer layer and into the log layer for the case of unstable stratification (figure 4c,f).
(a–c) Interface topography between the fluid and the solid at the final time of the simulation. Channels in the ice are highlighted in brown while keels are highlighted in green. (d–f) Spanwise variations of interface topography 
$\xi '=\tilde {\xi }-\bar {\xi }$ normalized by the viscous length scale 
$\delta _{\nu }=1/Re_*$ and as functions of 
$t_*$. The top, middle and bottom rows show the results for stable, neutral and unstable stratification, respectively.
Figure 4(d–f) shows that the viability of channels and keels increases with decreasing stratification: channels and keels are short lived with stable stratification but long lived with unstable stratification. For stable stratification (figure 4d), the separation of scales between the topography lifetime (about 10 friction time units) and the diffusion time scale across the solid layer (about 100 friction time units) suggests that the interface evolution is purely driven by the flow dynamics. For neutral stratification, figure 4(e) shows that channels and keels can drift, merge, split, decay and spontaneously appear over time scales of tens to hundreds of friction time units, highlighting a possible interplay between interface evolution and the fixed-temperature condition at the top solid boundary. For unstable stratification (figure 4f), the channels and keels become time invariant and their amplitudes saturate because of the top solid boundary condition, which plays a key role in the interface evolution as discussed in the next sections.
3.3. Coupled dynamics of the fluid and solid phases
The emergence of channels and keels can be the result of either (i) a passive response of the interface to uneven melting patterns driven by the turbulent flow, or (ii) a fully coupled interplay between fluid turbulence, interface topography and temperature in the solid. Here, we investigate the relevance of regimes (i) and (ii) for each of our simulations by looking at both the flow dynamics and the temperature field in the solid.
Figure 5 shows the 
$x$-averaged vertical heat fluxes in both the fluid and the solid for stable (top panel), neutral (middle panel) and unstable stratification (bottom panel) at a late time when a statistical steady state is reached. In the lower figure of each panel, we show the heat map of the vertical heat flux 
$\tilde {q}=\widetilde {wT}-\partial _z\tilde {T}$ in the fluid and of the conductive heat flux 
$-\partial _z\tilde {T}$ in the solid (note that we subtract 
$Nu^{III}$ in order to highlight fluctuations and that 
$-\partial _z\tilde {T}$ provides a more accurate measure than 
$\tilde {q}$ for the heat flux in the solid; cf. appendix B). In all three simulations, the spanwise fluctuations of the vertical heat flux are one order of magnitude (or more) larger in the fluid bulk where convection is active than in the solid where there is no movement. This suggests that the fluid flow has the capacity to induce transient melting or freezing hot spots and generate topography on short time scales, whereas temperature diffusion in the solid is moderate and primarily passive since it acts on long diffusive time scales of the order of 100 friction time units. The upper figure of each panel shows the melt rate 
$\tilde {\dot {m}}$ as a function of 
$y$ (dotted blue line), as well as the vertical conductive heat fluxes in the fluid (solid red line) and solid (dashed green line) right below and above the interface, as illustrated in the top right panel. As expected, the melt rate is positive when the heat flux coming from the fluid exceeds the heat flux going through the solid, i.e. when the red solid line is above the green dashed line. In figure 5(c), for instance, melting occurs on the left-hand side of each of the two channels (
$y\approx 2.5$ and 
$y \approx 5.6$) because the local, instantaneous heat flux coming from the fluid (red solid line) exceeds the heat flux going through the solid (green dashed line). Note, however, that melting and freezing events balance each other on long time scales in the simulation with unstable stratification such that, in agreement with figure 4(f), the channels do not migrate in the spanwise direction. The spanwise fluctuations of 
$-\partial _z\tilde {T}$ are moderately (respectively slightly) larger in the fluid than in the solid for stable (respectively neutral) stratifications. This suggests that the fluid dynamics controls the topography for stable and neutral stratification. In fact, we can observe that most local maxima of 
$-\partial _z\tilde {T}$ in the fluid near the boundary shown by the red solid lines in the line plots of figure 5 correspond to large upwelling events in the fluid bulk (see e.g. 
$y\approx 4$ for the top panel and 
$y\approx 0$ for the middle panel). However, heat fluxes in the solid still play a key role since there would be no freezing for a solid with a uniform temperature field. Heat flux fluctuations increase in the solid phase as well as in the fluid phase with decreasing 
$Ri_*$ (i.e. from top to bottom row): enhancement of heat flux fluctuations are the result of buoyancy effects in the fluid and of larger interface deformation in the solid. All heat maps show that 
$-\partial _z\tilde {T}-Nu^{III}$ in the solid is generally positive (green) over a channel and negative (pink) over a keel. In the case of stable and neutral stratification, 
$-\partial _z\tilde {T}-Nu^{III}$ decays from large fluctuation values near the interface to almost 0 near the top boundary, suggesting no significant influence of the top boundary on the temperature field in the solid (i.e. increasing the ice thickness while adjusting the top temperature to conserve the heat flux would not change the results). On the other hand, in the case of unstable stratification (figure 5c), 
$-\partial _z\tilde {T}$ fluctuations remain large near the top boundary, suggesting that there is a backreaction from the fixed-temperature top-boundary condition, 
$T=-q^s/2$ at 
$z=1.5$, on the interface evolution. The backreaction from the top boundary for unstable stratification is obtained because the position of the channels and keels becomes rapidly stationary (in contrast to the stable and neutral cases), such that the temperature field in the solid has time to adjust diffusively and balance the growth of the channels and keels (cf. appendix D for more details on the temperature field in the solid). The steadiness of the fluid dynamics and interface topography for unstable stratification can be seen to result in overlapping heat fluxes in the top figure of the bottom panel.
Panels of vertical heat fluxes for (a) stable, (b) neutral and (c) unstable stratification at time 
$t_*\approx 500$ representative of the statistical steady states. The lower figure of each panel shows the 
$x$-averaged vertical heat fluxes 
$-\partial _z\tilde {T}-Nu^{III}$ in the solid (pink–green colourmap) and 
$\tilde {q}-Nu^{III}$ in the fluid (spectral colourmap), the 
$x$-averaged interface position 
$\tilde {\xi }$ (black solid line), and the velocity vectors 
$(\tilde {v},\tilde {w})$ in the fluid (black arrows) and the conductive heat flux 
$(-\partial _y\tilde {T},-\partial _z\tilde {T})$ in the solid (red arrows). The upper figure of each panel shows the vertical conductive heat flux 
$-\partial _z\tilde {T}-Nu^{III}$ at a distance 
$3\delta _{\nu }=0.02$ above (green dashed lines) and below (solid red lines) the interface, and the melt rate (blue dotted lines), as illustrated in the upper right panel where the grey shading highlights the thickness of the diffuse interface.
The emergence of streamwise channels and keels is consistent with the well-documented presence of near-wall streamwise streaks and vortices in stratified shear flows (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017; Zonta & Soldati Reference Zonta and Soldati2018), but their amplitude clearly varies with 
$Ri_*$. In the case of stable stratification (figure 5a), buoyancy effects inhibit the generation of large topographical features, such that channels and keels have small amplitudes and do not feed back onto the flow (which is further discussed in § 3.4). In the case of neutral stratification, buoyancy is turned off, such that the solid boundary deforms more and can affect the flow dynamics. Figure 5(b) shows that the heat flux in the fluid close to the boundary is usually larger where there are channels (e.g. 
$y\approx 1.5,2.5$) than where there are keels (e.g. 
$y=0.9$), suggesting a topographic influence on the flow. For the case of unstable stratification (figure 5c), two streamwise rolls aligned with the direction of the flow and filling the entire depth dominate the fluid dynamics. These flow features are reminiscent of Rayleigh–Bénard convection rolls as observed in channel flow simulations with unstable stratification (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017), which here appear locked within the interface deformation pattern (figure 5c).
The interplay between the solid boundary and the flow dynamics for the case of unstable stratification is further highlighted in figures 6 and 7. Figure 6 shows the Hövmoller diagram of the heat flux in the middle of the fluid (figure 6a) and at 3 wall units below (figure 6b) and above the interface (figure 6c). Two mid-depth streamwise rolls whose positions are locked are evident from 
$t_*\approx 370$ onward in figure 6(a), which is about the same time as when the two channels and keels become large in figure 4(f). Similar rolls can be inferred for 
$t_*<370$ but are weaker and meander. The two strong rolls for 
$t_*\geq 370$ are locked with the topography (cf. figure 4f) and support large conductive heat fluxes with similar patterns right below the interface (cf. figure 6b), which further demonstrates that global modes control the interface dynamics for an unstable stratification. The conductive heat fluxes coming from the fluid are yet eventually balanced by the conductive heat fluxes through the solid (figure 6c), which adjust diffusively as the interface deforms, such that there is no net melting (figure 6d) beyond the initial transient of stage III. The decrease of conductive heat flux below the interface between 
$t_*^{Ic}$ and 
$t_*^{II}$ in figure 6(b) is due to the heat flux imbalance at the beginning of stage II (cf. appendix A) and has no incidence on the subsequent melting dynamics. Note that the temporal resolution is relatively coarse in figure 6(a–d) because these figures required the full three-dimensional data sets (due to interpolations in 
$z$ for figure 6b,c), which we saved at a low frequency to minimize disk usage. Also, the patchiness of the melt rate in figure 6(d) is due to the interaction of the ice boundary with individual turbulent bursts, including small-scale eddies, which drive intermittent melting and freezing events but have no long-term effect on topography generation.
(a) Heat flux anomaly 
$q'=(wT)'-T_z'$ in the fluid at 
$z=0.5$ and conductive heat flux anomaly 
$-T_z'$ (b) in the fluid at 
$\chi =-0.02$ (3 wall units) below the interface and (c) in the solid at 
$\chi =+0.02$ above the interface as functions of 
$(t_*,y)$ for the simulation with unstable stratification. The two vertical dashed lines highlight the times at which we turn on volume penalization and melting. (d) Melt rate 
$\tilde {\dot {m}}$ averaged in the streamwise direction as a function of 
$(t_*,y)$.
Spanwise spectrum of the heat flux 
$\tilde {q}$ in the fluid at 
$z=0.5$ (solid lines) and at 
$\chi =-0.04$ below the interface (dashed lines). Black and red colours highlight the spectra obtained in stage I and III, respectively (with 100 friction time units averaging). The blue dotted line shows the spectrum of 
$\tilde {\xi }$ averaged over the last 100 friction time units of stage III.
Figure 7 shows the spanwise spectrum of 
$\tilde {q}$ (left axis) in the middle of the fluid (solid lines) and near the solid boundary (dashed lines) at the end of stages I (red) and III (black) for the unstable case. In all cases, the spectrum peaks at wavenumber 
$L_yk_y/(2{\rm \pi}) =2$, which is close to the critical wavenumber 
$L_yk_c/(2{\rm \pi}) = 3.11$ of convection instability (
$\lambda _c=2.016$), suggesting that Rayleigh–Bénard convection is already active in stage I. However, with melting turned on, the peak is significantly amplified, especially for the spectrum near the boundary (dashed lines), which also shows amplification of higher harmonics, consistent with the spectrum of the interface itself (dotted blue line; right axis). These results suggest that Rayleigh–Bénard rolls are energized more than any other fluid features once melting is turned on, because they best couple with the interface topography evolution as a result of melting and freezing. Note that the spectra of the heat flux near the boundary and of the interface have a sawtooth-like pattern due to the non-sinusoidal shape of the interface and numerical confinement in the spanwise direction.
We next show in figure 8(a–f) (top row) the vertical profiles of the mean streamwise velocity 
$U$ and of the root-mean-square (r.m.s.) vertical velocity, 
$w_{rms}$, where here the mean and r.m.s. are defined using a horizontal and temporal average. The results are shown for all three stages and for stable (two leftmost columns), neutral (two middle columns) and unstable (two rightmost columns) stratification. The vertical profiles are never symmetric with respect to the half-fluid depth position, which is 
$z=0.5$ in stages I and II. This is because our velocity boundary conditions across the fluid layer are different, i.e. no-slip on the top boundary (which can be the fluid–solid interface) and free-slip on the bottom. For stable and neutral conditions, there is a strong overlap of all curves, suggesting that the topography does not influence the mean profiles, while for unstable stratification, there is a small deviation of the stage III profiles (solid lines). Note that a log–log plot of the mean velocity and temperature profiles near 
$z=1$, which we show in figure 14 in appendix A, clearly shows that the flows follow the law of the wall on the top no-slip fluid boundary in stages I and II.
Vertical profiles of the streamwise velocity 
$U$ and vertical velocity r.m.s. 
$w_{rms}$ averaged in 
$x$ and 
$t_*$ (over at least 50 friction time units) for stable (two leftmost columns), neutral (two middle columns) and unstable (two rightmost columns) stratification. The top panels show the full horizontal average of 
$U$ and 
$w_{rms}$, i.e. when they are also averaged in 
$y$, and the mean interface position 
$\bar {\xi }$ at the end of stage III (dotted lines). The bottom panels show the full horizontal averages of 
$U$ and 
$w_{rms}$ as well (dotted lines), but also 
$U$ and 
$w_{rms}$ under the largest keel (downward triangles) and channel (upward triangles), i.e. they are not averaged in 
$y$, as functions of 
$1-\chi$, where 
$\chi$ is the distance from the interface. In (
$l$), the green (respectively blue) thick and thin dashed lines show the vertical profiles shifted away from the largest keel (respectively channel) by 0.25 and 0.5 units in the 
$+y$ direction.
We investigate in figure 8(g–l) whether the mean profiles vary in the spanwise direction in such a way that they correlate with the 
$x$-averaged interface topography. To do so we consider the mean profiles under the largest keel (downward triangles) – i.e. at 
$y=y_k$, where 
$y_k$ is the minimum in 
$y$ of 
$\tilde {\xi }$ at each time step – separately from the mean profiles under the largest channel (upward triangle) – i.e. at 
$y=y_c$, where 
$y_c$ is the minimum in 
$y$ of 
$\tilde {\xi }$ at each time step – where mean now denotes streamwise and temporal averaging. or stable stratification, there is no difference between the profiles under keels and channels. However, for neutral and unstable stratification, the two profiles depart in such a way that the streamwise velocity is larger under channels than under keels, and the vertical velocity r.m.s. (defined using a temporal and 
$(x,y)$ average for each 
$\chi$) is larger under keels than under channels. These results indicate a noticeable influence of the topography on the flow. For the case of neutral stratification, the separation of the profiles is maximum for 
$\chi <0.3$ and then vanishes, suggesting a local influence of the topography on the flow dynamics, while for unstable stratification, the effect of the topography is felt throughout the entire depth due to coupling with the Rayleigh–Bénard rolls. It may be noted that the profiles of 
$w_{rms}$ under the largest keels and channels are larger than the plane-average profile shown by the dotted line in figure 8(l). This is expected because Rayleigh–Bénard convection promotes both localized intense upwellings and intense downwellings under channels and keels. In fact, away from the main channel and keel the profiles decrease rapidly, as can be seen from the blue and green lines.
In order to gain further insight into the statistics of the flow interacting with the melting boundary, we show in figure 9 the probability density functions (p.d.f.s) of the streamwise velocity (a,d,g), the vertical velocity (b,e,h) and the temperature gradient (c,f,i), for stable (a–c), neutral (d–f) and unstable (g–i) stratification. For the velocities, the p.d.f.s are shown both in the middle of the fluid, at 
$z=0.5$, and near the boundary, at 
$z=\xi -0.04$ (i.e. 6 wall units into the fluid). We find little difference between stages I (dashed lines) and III (solid lines) for the streamwise and vertical velocities, suggesting limited influence of the topography on the overall flow morphology, although the streamwise velocity in figure 9(g) has a negative tail with higher probability density in stage III than in stage I. The temperature gradient at the fluid–solid interface (figure 9c,f,i) does vary noticeably between stage I and stage III. However, this difference is due to the phase-field method rather than to a fundamental change in flow morphology since the p.d.f.s in stages III and II (not shown) show significant overlap.
Probability density functions of the streamwise velocity fluctuations 
$u-\bar{{\bar {u}}}$ (leftmost panels), vertical velocity fluctuations 
$w-\bar{{\bar {w}}}$ (middle panels) and temperature gradient fluctuations 
$\partial _zT-\overline {\overline {\partial _zT}}$ (rightmost panels), for stable (top row), neutral (middle row) and unstable (bottom row) density stratification, with the double overbar denoting the horizontal and temporal average. The p.d.f.s of 
$u-\bar{{\bar {u}}}$ and 
$w-\bar{{\bar {w}}}$ are shown at 
$z=0.5$ and 
$z=\xi -0.04$ (6 wall units below the interface), while the p.d.f.s of 
$\partial _zT-\overline {\overline {\partial _zT}}$ are shown at the interface 
$z=\xi$. We use dashed lines and solid lines to highlight the results obtained in stage I and stage III, respectively. Blue and red represent the p.d.f.s of 
$u-\bar{{\bar {u}}}$ and 
$w-\bar{{\bar {w}}}$ conditioned on negative and positive local temperature fluctuations, in order to separate the p.d.f.s of warm fluid patches from those of cold fluid patches.
While phase changes and the emergence of topographical features have little effect on the p.d.f.s, most p.d.f.s display flow-driven left–right asymmetries, which are worth highlighting. Most importantly, the p.d.f. of the temperature gradient at the fluid–solid interface has a rapidly decaying positive tail and a slowly decaying negative tail. This asymmetry is obtained in all stages and hence is a feature of the flow rather than a consequence of topography generation, and suggests that the phase-change dynamics should be itself asymmetric (which we show in § 3.5). The p.d.f.s of the streamwise velocity near the boundary are also asymmetric, featuring a slowly decaying positive tail and a rapidly decaying negative tail. We have further separated the p.d.f.s of the velocities in figure 9 based on the sign of the local temperature anomaly (compared with the plane and temporal mean). Blue curves denote p.d.f.s obtained for a negative temperature anomaly, i.e. representative of fluid patches influenced by the cold top boundary, while red curves denote p.d.f.s obtained for a positive temperature anomaly, i.e. representative of fluid patches influenced by the warm bottom boundary. The cold temperature p.d.f.s are shifted to the left of the warm temperature p.d.f.s for the streamwise velocity (left column), which suggests that negative streamwise velocity is more often associated with cold fluid coming from the top boundary. Also, for the vertical velocity at 
$z=\xi -0.04$ (right panels of the middle column), the positive tail of the warm temperature p.d.f.s (red) is larger than the negative tail of the cold temperature p.d.f.s (blue), suggesting more extreme warm upwelling events than cold downwelling events just outside of the viscous sublayer. These results demonstrate that the near-wall flow dynamics have multiple asymmetries, which may be related to the asymmetry in the temperature gradient at the boundary.
3.4. Reversing the stratification
In the case of stable stratification, the cool melt fluid (at freezing temperature) is more buoyant than the warmer surrounding fluid, such that it rises and accumulates in the middle of the channels. Thus, channels and keels are limited to small amplitudes in figure 5(a) because of a negative feedback between channel generation and the melt dynamics, which halts channel growth. In the case of unstable stratification the opposite is true, i.e. the cool melt fluid is denser than the surrounding fluid and, hence, is evacuated from channels. This yields a positive feedback between channel generation and melt dynamics, which results in the large-scale topographical features seen in figure 5(c) (which saturate over time due to thermal adjustment in the solid).
The hypothesis of the cool melt fluid pooling in channels and inhibiting their growth for stable stratification is difficult to verify with the results discussed previously because of the small interface deformation obtained for 
$Ri_*=40$. Therefore, we have run a fourth simulation starting from the final time of the simulation with unstable stratification (and large interface deformation), but with an increasing Richardson number such that the fluid becomes stably stratified and interacts (transiently) with the initially large topographical features. We impose the stable stratification through several intermediate steps so that the flow does not relax to a laminar state. Specifically, we use
\begin{equation} Ri_*(t_*) = -40 [1 - f(t_*,5)] + 20 [2 + f(t_*,12) + f(t_*,19)], \end{equation}
with 
$f(t_*,\tau _*)=\tanh ( t_*-t_*^{III} - \tau _*)$, such that 
$Ri_*$ starts from 
${\approx }{-}40$ at 
$t_*=t_*^{III}$ and reaches 
${\approx }40$ for 
$t_*>t_*^{III}+19$. We show the results of this run in figure 10, where blue/red highlight 
$x$-averaged temperature values in the solid/fluid phase, while arrows denote 
$x$-averaged velocity vectors 
$(\tilde {v},\tilde {w})$ in the 
$(y,z)$ plane. Figure 10(a) shows the results at time 
$t_*=t_*^{III}+1$ (
$t_*^{III}=526$), i.e. when 
$Ri_*\approx -40$. The stratification is unstable such that the flow features strong upwelling of warm fluid below the channels and strong downwelling of cool fluid along and under the keels, akin to Rayleigh–Bénard rolls locked into the deformed interface pattern. A pair of counter-rotating streamwise rolls is clearly visible below each of the two channels. These rolls persist until 
$Ri_*\approx 0$, which is in agreement with recent simulations of mixed convection that have shown that streamwise rolls extending throughout the full depth of a channel (without phase changes) are obtained for a wide range of negative Richardson numbers (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). Figure 10(b,c) show the results at times 
$t_*=t_*^{III}+13$ and 
$t_*=t_*^{III}+21$, i.e. when 
$Ri_*\approx 0.7$ and 
$Ri_*\approx 39$, respectively. At these times, the stratification is stable and the cool melt fluid produced at the keels converges toward the channels’ centreline. The Rayleigh–Bénard large-scale rolls have vanished and are replaced with weaker vortices of finite vertical extent, which are most vigorous close to the interface where they are driven by the (positive) buoyancy anomaly of the melt fluid at the tip of the keels. The heat flux through the fluid goes down and freezing occurs everywhere such that the solid front advances into the fluid. Importantly, freezing is faster in the channels because of the convergence of the buoyant cool melt fluid and higher conductive heat fluxes in the solid (as the solid is thin above channels), which leads to rapid refreezing of the initial channels. While 
$Ri_*>0$ increases, the properties of the boundary-attached vortices (e.g. vertical extent and intensity) are the result of a complex interplay between the amplitude of the interface topography and of the stratification strength. The increasing stratification drives an increasing positive buoyancy anomaly of the melt fluid but also increasingly damps global modes (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011), while the decreasing topography amplitude is expected to result in flattening and weakening vortices. Ultimately, the topography disappears and the vortices weaken significantly.
Temperature field averaged in the 
$x$ direction in the fluid (red colourmap) and in the solid (blue colourmap) at times (a) 
$t_*=527$, (b) 
$t_*=539$ and (c) 
$t_*=547$ for the restratifying simulation. The black arrows show the velocity field 
$(\tilde {v},\tilde {w})$. The thick solid line shows the 
$x$-averaged interface position 
$\tilde {\xi }$.
3.5. Melting
In this section we discuss the evolution of the mean interface position 
$\bar {\xi }$ with time and the statistics of melting 
$\dot {m}=\partial _t\xi$ at the statistical steady state. We first show the evolution of 
$\bar {\xi }$ in DNS as a function of time in figure 11(a) for all three simulations (solid lines). For a solid with spatially uniform temperature equal to the melting temperature 
$T_m$, we would expect a faster increase of 
$\bar {\xi }$ with time for an unstable than for a stable stratification, since the heat flux in the fluid is larger when buoyancy forces are destabilizing rather than restoring. However, as indicated in § 2, we have imposed a conductive heat flux in the solid (
$q^s$) slightly less than the mean heat flux through the fluid (
$Nu^I$) at the beginning of stage III, i.e. when we turn on melting, such that the leading-order melt rate is not controlled by the Nusselt number of the fluid-only simulations but by the difference 
$Nu^I-q^s$. This difference is 0.19, 0.1 and 0.05 for stable, neutral and unstable stratification (cf. table 1), such that the initial increase of 
$\bar {\xi }$ is faster for stable than for unstable stratification. Figure 11 shows that 
$\bar {\xi }$ saturates over time. This happens because the conductive heat flux in the ice increases as 
$\bar {\xi }$ increases (in a plane-averaged sense), since the solid becomes thinner, such that it eventually balances the heat coming from the fluid.
(a) Time history of the mean interface position for stable, neutral and unstable stratification (thick solid lines). The dotted (respectively dashed) lines show the position of the interface predicted by the reduced model (3.4) with 
$Nu^I$ (respectively 
$Nu^{III}$) as the heat flux in the fluid. The dash–dot line shows the result of the higher-order model (D1), (D10) for unstable stratification. (b) Time history of the heat flux relative to 
$Nu^I$, i.e. 
$(Nu-Nu^{I})/Nu^I$ (solid lines). We use a rolling mean over a window of 
$\varDelta _*=20$ friction time units in order to remove some of the most rapid large-amplitude oscillations. The horizontal dashed lines show 
$(Nu^{III}-Nu^{I})/Nu^I$.
Under the assumption of small interface deformations, it is possible to predict the evolution of the mean interface position over time using a reduced model. As a first approximation, we consider that the topography has no effect on the temperature in the solid, i.e. we assume that the heat flux through the solid is simply equal to the temperature difference between the interface and the top boundary divided by the mean solid thickness 
$h$ (cf. details in appendix D, which hinges on an assumption of a quasi-steady state). Then, the evolution equation for 
$\bar {\xi }$ becomes
\begin{equation} \frac{\textrm{d}\bar{\xi}}{\textrm{d} t} = q^f - \frac{h_0q^s}{h}, \end{equation}
where 
$h_0=1/2$ is the initial ice thickness, 
$q^f$ is the heat flux in the fluid and we recall that 
$\bar {\xi }=\xi _0=1$ at 
$t=t_*^{II}$. For simplicity, we take 
$q^f$ to be a constant diagnosed from the simulations. The results of (3.4) for 
$q^f=Nu^I$, i.e. obtained when setting 
$q^f$ to the average heat flux before melting is turned on, are shown by the dotted lines in figure 11(a). The overlap between the reduced model and DNS results at early times is good for unstable stratification (as expected) but is poor for stable and neutral stratification. The disagreement with 
$q^f=Nu^I$ at early times arises because the temperature in the solid is slightly above 0 because of volume penalization, such that there is some artificially large melting at the beginning of stage III (cf. appendix A). At later times, the DNS results and the model results shown by the dotted lines diverge because the heat flux 
$Nu$ increases rapidly once melting is turned on, as can be seen in figure 11(b). For unstable stratification, the agreement with 
$q^f=Nu^I$ is relatively good until 
$t_*\approx t_*^{II}+50$ (cf. red dotted line), i.e. right until the Rayleigh–Bénard rolls are energized and a large-scale topography emerges (cf. § 3.3).
In order to account for the increase in heat flux through the fluid enabled by melting and the generation of topography, we show with dashed lines in figure 11(a) the result of (3.4) with 
$q^f=Nu^{III}$, which is the heat flux at the statistical steady state with melting turned on. For stable and neutral stratification, there is good agreement between the model results and the mean interface position at late times. For unstable stratification, however, (3.4) with 
$q^f=Nu^{III}$ overestimates the final value of 
$\bar {\xi }-1$ by a factor of two (approximately), suggesting that topography plays a non-negligible role on the heat flux in the solid. We show in figure 11(a) a prediction of 
$\bar {\xi }$ for unstable stratification obtained using a more accurate higher-order model (red dash–dot line), which takes into account interface deformation (cf. appendix D). The higher-order prediction overlaps well with the DNS results at late times, demonstrating that melting and the generation of topography changes the heat flux through both the fluid and the solid. The topography makes the solid more efficient at evacuating heat because the anomalous (increased) heating obtained above the channels (i.e. where the solid is thin) exceeds in absolute value the anomalous (reduced) heating obtained above the keels (i.e. where the solid is thick), which is a nonlinear effect in the topography amplitude obtained for any topography with top-down symmetry (e.g. a sinusoid). The higher-order model takes into account this nonlinear effect in topography amplitude and predicts a steady-state solid layer thickness larger than that predicted by the low-order model without topography for the same forcing heat flux 
$q^f$.
We finally show in figure 12 the p.d.f.s of interface deformation and melt rate at the statistical steady state, i.e. past the initial transient during which topographical features emerge and the solid melts on average. At the statistical steady state, the p.d.f.s of interface position become roughly time invariant. The mean interface position reaches a plateau (cf. figure 11a) because the mean amount of freezing balances the mean amount of melting at every time step, and the standard deviation, or topography amplitude, saturates (as can be seen in figure 4(f) for the unstable case). For stable stratification, the steady-state p.d.f. of interface deformation is almost symmetric with respect to the mean and appears approximately Gaussian (figure 12a). This suggests that channels and keels are symmetric to each other with respect to the mean interface position for the stable case. For neutral stratification, a small asymmetry develops, i.e. the median shifts toward small positive deformations (channels) and the tail of extreme negative deformations (keels) increases slightly. The same asymmetry is amplified for unstable stratification, with a narrow peak appearing to the right of the 0 mean and the negative tail increasing further. In other words, as the stratification becomes unstable, patterns grow in size and the width-to-height ratio of channels increases (broad and flat) while the width-to-height ratio of keels decreases (narrow and deep). The asymmetry in the p.d.f.s of interface topography is consistent with the observation from figure 4 that channels are typically flatter and more widespread than keels.
Probability density functions of (a) the interface position 
$\xi$ and (b) the melt rate 
$\dot {m}$, for neutral, stable and unstable stratification. Double overlines denote temporal (100 friction time units) and spatial averaging. In (b) we show the p.d.f.s for the full melt rate (solid lines), the melt rate along channels (upper triangles) and the melt rate along keels (lower triangles).
The p.d.f.s for the melt rate 
$\dot {m}$ are shown in figure 12(b). The temporal and spatial average, 
$\bar {\bar {\dot {m}}}$, which is subtracted from the p.d.f., is close to 0 in all cases, since the mean amount of melting is balanced by the mean amount of freezing at the statistical steady state. The p.d.f.s of melt rate are asymmetric, i.e. similar to the p.d.f.s of heat flux at the top of fluid-only channel simulations (see dashed lines in figures 9c,f,i). The median is shifted to the left of the mean, i.e. toward negative values representative of freezing events, and the positive tail is enhanced compared with the negative tail. In other words, the interface is typically freezing slowly (
$\dot {m}<0$), but occasionally melts rapidly (
$\dot {m}>0$). We remark that the asymmetry of the melt rate p.d.f.s is not due to the asymmetry of the interface p.d.f.s since the p.d.f.s of melt rates inside channels (upper triangle) and along keels (lower triangles) are similar, but is instead a generic feature of melting by a turbulent fluid. Indeed, while the turbulent flow can drive rapid melting independently of what happens in the solid, freezing necessarily involves slow diffusive processes in the solid. Additionally, the near-wall dynamics, which features coherent structures such as streamwise streaks and vortices, is itself asymmetric, as can be seen from the p.d.f.s of the temperature gradient in figure 9. Thus, it is not surprising that the melt rate p.d.f.s are asymmetric. While beyond our goal, it would be worthwhile in the future to try to identify flow features controlling the shape of the melt rate p.d.f.s.
4. Geophysical discussion
Due to the large computational costs of coupled fluid–solid simulations, all control parameters in this study were held fixed, i.e. we considered 
$Re_*=150$, 
$Pr=1$ and 
$St=1$, except for the friction Richardson numbers, which we varied in order to test the effect of density stratification on topography generation. Our simulation with a positive (respectively negative) Richardson number 
$Ri_*=40$ (respectively 
$Ri_*=-40$), i.e. with stable (respectively unstable) stratification, assumes a negative (respectively positive) thermal expansion coefficient. Cold freshwater has a negative (respectively positive) thermal expansion coefficient at low (respectively high) pressure (Thoma et al. Reference Thoma, Grosfeld, Smith and Mayer2010). Thus, our stable simulation is qualitatively similar to the flow of freshwater below an ice cover at low pressure (as is the case in an ice-covered lake), whereas our unstable simulation is qualitatively similar to the flow of freshwater below an ice cover at high pressure (as is the case in a deep subglacial lake). In the case where the solid is below the fluid, the stratification is reversed, i.e. the unstable simulation results are applicable to the flow of cold freshwater above ice at low pressure (as is the case in supraglacial rivers). Our stable simulation is also qualitatively similar to the flow of salt water under ice shelves. The melt water under ice shelves is cooler but also fresher than the ambient ocean water, such that it is positively buoyant. In fact, salinity and temperature can be combined, assuming that they have the same effective diffusivities, into a single variable known as thermal driving, which has a negative expansion coefficient (Jenkins Reference Jenkins2016).
In order to minimize resolution requirements and observe large topographical changes in a relatively small amount of time, we have considered a flow that is only moderately turbulent, weakly stratified and anomalously warm. If we assume that the working fluid is water, i.e. 
$\mathcal {L}=3\times 10^5$ J kg
$^{-1}$ and 
$c_p=4\times 10^3$ J kg
$^{-1}$ K
$^{-1}$, then 
$St=1$ implies that the bottom boundary is held at 
$75\,^{\circ }\mathrm {C}$. If we further assume that 
$H=10$ cm and 
$\nu =\kappa =10^{-6}$ m
$^2$ s
$^{-1}$, i.e. considering an anomalously high thermal diffusivity such that 
$Pr=1$, then our simulation runs correspond to about 10 hours in real time, the bulk velocity is 2 cm s
$^{-1}$ and the thermal expansion coefficient is 
$5\times 10^{-7}$ K
$^{-1}$ in absolute value, which is quite small for water. These calculations highlight that there is clearly a significant gap between the control parameters in our simulations and in nature. Our set of experiments is also limited to three different runs, such that we cannot offer a quantitative prediction of what would be observed in the environment. Nevertheless, Rayleigh–Bénard rolls have been observed for a wide range of Rayleigh numbers and Reynolds numbers in mixed-convection simulations with unstable stratification (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017), such that large-scale channels and keels driven by these rolls may be expected in most conditions with unstable stratification, for example in supraglacial rivers or deep subglacial lakes, especially when the external flow is weak. The transverse wavelength of channels and keels maintained by global rolls is expected to be of the same order as the fluid depth, as is the case in our unstable simulation (cf. figure 10(a), in which each wavelength accommodates two counter-rotating rolls with a diameter equal to the mean fluid depth). We note that our unstable simulation has low 
$Ra=4.5\times 10^5$ and 
$Re_*=150$ compared with state-of-the-art simulations of mixed convection (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020), but more importantly has a relatively small bulk Richardson number 
$Ri_b=RaPr/Re_b^2\approx 0.1$ (defined as positive for an unstable stratification) and large Monin–Obukhov length 
$L_{MO}\approx 1$ (normalized by 
$H$). Thus, channels and keels can be expected for a broad range of bulk Richardson numbers of unstably stratified shear flows, i.e. 
$Ri_b\geq 0.1$, and to grow in size as 
$Ri_b$ increases. In the limit 
$Ri_b \rightarrow \infty$, channels and keels may be expected to eventually disappear. Indeed, as 
$Ri_b \rightarrow \infty$, the mean flow vanishes and buoyancy effects dominate, such that three-dimensional domes and cusps should emerge in place of streamwise features (Rabbanipour Esfahani et al. Reference Rabbanipour Esfahani, Hirata, Berti and Calzavarini2018). Interestingly, large channels have also been observed at the base of ice shelves. However, these channels are unlikely to originate from Rayleigh–Bénard rolls but rather from transverse perturbations of, for example, the subglacial discharge or ice thickness at the grounding line (Dallaston, Hewitt & Wells Reference Dallaston, Hewitt and Wells2015), since the stratification is in this case stable.
For stable and neutral stratifications, we also observed channels and keels. Channels and keels with stable stratification (
$Ri_*\geq 0$) are carved by boundary-attached momentum streaks rather than by global modes, however, such that their transverse wavelength is shorter than for the case of unstable stratification (although here the difference is weak given the small 
$Re_*$), and their amplitude is either equal to or smaller than the viscous sublayer thickness, i.e. small. The shape and size of the small channels and keels obtained for 
$Ri_*\geq 0$ are in stark contrast with the three-dimensional scallops, which have been observed in stably stratified polar oceans and neutral laboratory experiments. Scallops observed in the field and investigated in laboratory experiments have amplitudes of the order of a few centimetres and wavelengths of the order of a few tens of centimetres, i.e. they are tall and wide features compared with the viscous sublayer thickness, which is typically smaller than 1 mm in nature (Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019). Previous experimental and theoretical works have found that the friction Reynolds number based on the scallop wavelength 
$\lambda$ usually satisfies 
$Re_*^{\lambda }=Re_*\lambda /H \geq O(1000\text {--}10\,000)$ and have always reported a scallop wavelength smaller than the fluid depth, i.e. 
$\lambda < H$ (Blumberg & Curl Reference Blumberg and Curl1974; Thomas Reference Thomas1979; Claudin et al. Reference Claudin, Durán and Andreotti2017). Considering the upper limit 
$\lambda = H$ means that scallops are predicted to emerge for 
$Re_*\geq O(1000\text {--}10\,000)$ in water, which is at least one order of magnitude higher than what we selected for our study and difficult to achieve numerically. Note, however, that the minimum 
$Re_*$ leading to scallop formation may be different for a fluid with control parameters 
$Pr=1$ and 
$St=1$ (as is the case in this work) instead of 
$Pr\approx 10$ and 
$St\approx 75$, as is the case for water, and may also vary with the stratification strength.
5. Concluding remarks
We have shown that streamwise channels and keels spontaneously emerge as the dominant topographical features of a fluid–solid boundary when the flow is pressure-driven, turbulent and thermally stratified with 
$Re_*=150$, 
$Pr=1$ and 
$St=1$. We have investigated the effect of the background density stratification and found that the amplitude of the channels and keels increases with decreasing stratification. For unstable stratification (
$Ri_*=-40$), the channels and keels couple strongly with Rayleigh–Bénard rolls, which are energised and locked within the interface deformation pattern. For neutral stratification, a similar correlation is obtained between the flow dynamics and the interface deformation pattern. However, the full-depth rolls are replaced with smaller and weaker boundary-attached momentum streaks, which do not provide a clear locking mechanism, i.e. the topography drifts. For neutral (
$Ri_*=0$) and stable stratification (
$Ri_*=40$), the channels and keels saturate either because of the absence of a positive feedback between topography and momentum streaks or because stabilizing buoyancy forces inhibit channel growth. For unstable stratification, the saturation is due to the fact that we impose the temperature at the top boundary. With an imposed heat flux at the top, the entire solid would melt rapidly and entirely provided that the stratification is unstable (not shown), which means that the choice of boundary conditions at the top of the solid can be critical. Note that the growth of the fluid layer for unstable stratification is due to the positive feedback that melting has on the effective Rayleigh number of the convective fluid. As the solid melts, the effective Rayleigh number increases, leading to further melting, which is stopped only if diffusion in the solid can eventually balance the increasing heat flux in the fluid.
The analysis of the melt rate statistics indicates that there is an asymmetry in melting and freezing, which may be related to the different melting/freezing dynamics (freezing relying primarily on slow diffusive processes in the solid) but also asymmetries in the flow statistics. Specifically, melting is highly localized and intense while freezing is widespread but weak. While beyond the scope of this study, it would be useful to identify whether coherent features of the near-wall turbulent flow, such as streamwise streaks and vortices, correlate preferentially with either melting or freezing events.
The drag coefficient changes significantly depending on the type of stratification but is only weakly affected by the generation of topographical features, which is not unexpected in our case since streamwise channels and keels are smooth in the direction of the flow. Capturing three-dimensional topographical features, such as scallops, which do affect momentum and heat transfers (Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019), in coupled fluid–solid simulations would be a major achievement, which could complement fluid-only simulations at planar ice boundaries (Gayen, Griffiths & Kerr Reference Gayen, Griffiths and Kerr2016; Keitzl, Mellado & Notz Reference Keitzl, Mellado and Notz2016a,Reference Keitzl, Mellado and Notzb; Mondal et al. Reference Mondal, Gayen, Griffiths and Kerr2019; Vreugdenhil & Taylor Reference Vreugdenhil and Taylor2019). However, as discussed in § 4, scallops may require much higher Reynolds numbers to form than 
$Re_*=150$. In fact, the minimum 
$Re_*$ for scallops could be too high for a phase-field method on most supercomputers. The cheapest test for evaluating the minimum Reynolds number leading to scallops would be to start the simulations with a longitudinally wavy boundary and investigate the initial evolution. The runtime would be reduced to a minimum. However, a high resolution (higher than say 
$1024^3$ with a spectral code) would still be required. It is noteworthy that simulations of a pure fluid at a fixed wavy boundary would already be useful in helping to verify or refine the most recent theoretical predictions on scallop formation and saturation (Claudin et al. Reference Claudin, Durán and Andreotti2017) and estimate the effect of the stratification strength, which has not yet been considered. We remark that capturing scallops in a water environment would require not only higher 
$Re_*$ but also higher 
$St$ and 
$Pr$, which would both incur significant computational overheads. Higher 
$Pr$ results in thinner thermal boundary layers, which could impact the near-wall dynamics and, for example, the asymmetry between melting and freezing. Higher 
$St$ results in slower melt rates, which could significantly change how interface patterns couple with transient flow features. In the case of unstable stratification, we might still expect that Rayleigh–Bénard rolls couple with the interface deformation pattern for high 
$St$, since they are relatively stationary flow features at least in the strong shear regime (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). For neutral stratification, however, the interface evolves over time scales similar to those of the flow dynamics for 
$St=1$ (figure 4e), such that increasing 
$St$ might significantly decrease the sensitivity of the interface topography to fluid anomalies. Finally, freshening effects, which are critical to ice–ocean interactions, would require adding slowly diffusing salt to the simulations, which constitutes yet another significant challenge for multi-phase DNS.
From a fundamental physics viewpoint, it would be interesting to investigate in detail how topographical features are modified when phase changes are driven by dissolution rather than by melting. The fluid–solid boundary conditions (Stefan condition) and scalar diffusivities are different between dissolution and melting experiments. However, similar longitudinal and rippled patterns have been observed in both cases (e.g. Allen (Reference Allen1971) for dissolution). It would be also worthwhile to explore the effect of phase changes and topographical features on the onset of global modes and the large-scale organization of mixed-convection flows, which are of interest to many fields of physics and engineering (Kelly Reference Kelly1994; Pabiou, Mergui & Bénard Reference Pabiou, Mergui and Bénard2005; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020). Finally, it would be useful to investigate potential analogies between ice patterns due to melting and freezing and the formation of sand ripples and dunes, which have been and continue to be extensively studied (e.g. Charru, Andreotti & Claudin Reference Charru, Andreotti and Claudin2013; Courrech du Pont, Narteau & Gao Reference Courrech du Pont, Narteau and Gao2014).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2020.1064.
Acknowledgements
We gratefully acknowledge fruitful discussions with members of the ISOBL group at the British Antarctic Survey and the University of Cambridge. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement 793450. We acknowledge PRACE for awarding us access to Marconi at CINECA, Italy.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Phase-field method
The phase-field method transforms the discontinuous two-phase two-domain problem into a continuous two-phase one-domain problem, which can be solved numerically using a pre-existing fluid code. In order to reproduce the original problem correctly, the resolution and parameters 
$A$, 
$B$, 
$C$ and 
$\varGamma$ of the phase-field equations (2.5) must satisfy several constraints (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020), which can be verified a posteriori by diagnosing the flow properties in the solid and fluid phases. Here, for our choice of resolution, we have 
$A=6/(5St)$, 
$B=(16/\delta ^2)\times 6/(5St)$ and 
$C=1$, with 
$\delta$ chosen such that it is equal to two times the local grid size at 
$z=1$ (initial interface position). Also, 
$\varGamma =(\delta /2.648228)^2$ and we require time steps to be always smaller than 
$\varGamma /2$.
We first assess the effect of the phase-field method and choice of parameters on the flow variables by showing in figure 13 the ratio of the kinetic energy averaged over the solid volume, 
$KE_s$, divided by the kinetic energy averaged over the fluid volume, 
$KE_f$. Figure 13 shows that 
$KE_s/KE_f<10^{-4}$ and that the fluctuations are within a factor of two of the mean, i.e. velocities in the fluid penetrate only very weakly into the solid and do not burst significantly.
Kinetic energy in the solid divided by kinetic energy in the fluid as a function of time 
$t_*$. Red, gold and blue colours denote results obtained for unstable, neutral and stable stratification, respectively.
We now further comment on the resolution requirements and our choice for the grid size and time step. We recall that 
$\delta$ is the thickness of the diffuse phase-field interface over which 
$\phi$ transitions from 1 in the fluid to 0 in the solid (see figure 1c). Here 
$\delta$ is an artificial length scale, such that it must be smaller than any physical length scale in the problem, while at the same time being larger than the grid size since it must be resolved numerically. For a boundary layer flow with 
$Pr=1$, the smallest length scale close to the interface is the viscous sublayer thickness, which is typically equal to a few times the viscous length scale 
$\delta _{\nu }=1/Re_*$. Here, we chose to have 
$\delta \approx \delta _{\nu }$ in all simulations, i.e. 
$\delta$ is equal to 1 wall unit 
${\rm \Delta} z^+=\delta$, such that the diffusive interface for the phase field is comprised within the viscous and thermal sublayers, as can be seen in figure 14. In order to resolve the diffusive interface, we selected a vertical Chebyshev basis with enough modes such that the collocation grid has a resolution 
$dz$ near the mean interface position equal to or less than 
$\delta _{\nu }/2$. The damping time scale 
$\varGamma$ in (2.5) is set to 
$\varGamma =(\delta /2.648228)^2$ in order to cancel first-order errors in the phase-field model (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020).
Wall-normalized velocity 
$U^+=\overline {Pr\,u/u_*}$ (circles, left axes) and temperature 
$T^+=\overline {u_*T/(PrNu)}$ (crosses, right axes) expressed in terms of wall units as functions of the wall coordinate 
$z^+=\overline {u_*(1-z)/Pr}$ for (a) stable, (b) neutral and (c) unstable stratification. Blue symbols denote results obtained for the simple channel flow configuration (i.e. for 
$t_*\in [t_*^{Ib},t_*^{Ic}]$) whereas orange symbols denote results obtained with volume penalization (i.e. for 
$t_*\in [t_*^{Ic},t_*^{II}]$). The plus symbols show the temperature profiles obtained with volume penalization but shifted to the right, i.e. into the fluid, by replacing 
$z^+$ with 
$z^++2.5$, 
$z^++1.16$ and 
$z^++0.78$ for stable, neutral and unstable stratification. The overbars denote horizontal averaging and time averaging over 20 friction time units. The vertical dotted lines show 
$z^+=-\delta ,0,\delta$, where 
$\delta$ is the diffuse interface thickness.
Figure 14(a–c) shows a semilog plot of the wall-normalized velocity 
$U^+$ (left axes) and wall-normalized temperature 
$T^+$ (right axes) as functions of wall units 
$z^+$ (
$z^+\geq 0$ denote positions in the fluid while 
$z^+<0$ denote positions in the solid) in stages I and II for stable, neutral and unstable stratification, respectively. In the viscous and thermal sublayers, which extend from 
$z^+=0$ to 
$z^+\approx 5$, we expect a linear scaling for both 
$U^+$ and 
$T^+$ with 
$z^+$, shown by the solid dashed lines. This linear scaling is perfectly satisfied by the DNS results in stage I (blue circles and blue crosses) as well as the DNS results in stage II (orange circles and crosses), except for 
$|z^+|<\delta$ (shown by the vertical solid lines), i.e. within the diffuse interface, which is expected since this is where the dynamics is artificially controlled by the phase-field equation. It may be noted that 
$T^+$ (orange crosses) is anomalously large for 
$z^+<\delta$ and in fact deviates from the true solution (blue crosses) slightly outside the diffuse interface. This discrepancy is due to the fact that the heat flux in the fluid is larger than the heat flux in the solid in stage II. The interface being fixed in stage II, the heat imbalance results in the heating of the solid, such that 
$T^+=0$ occurs at 
$z^+<0$ away from the fixed interface position 
$z^+=0$ (note that we use a symmetric logarithmic scale with a linear threshold at 
$|z^+|=0.1$). By shifting the temperature profile to the right such that 
$T^+=0$ is aligned with 
$z^+=0$ (red pluses), we recover a perfect linear scaling for the temperature within both the thermal fluid sublayer and the solid. Outside of the linear sublayer and the buffer layer, the mean vertical profiles exhibit a logarithmic behaviour.
Far from the top boundary, that is, for 
$z^+\approx 100$, 
$U^+$ shows a steeper scaling with 
$z^+$ for stable stratification than for neutral or unstable stratification. This is a consequence of buoyancy effects, which tend to decrease (respectively increase) stirring of the mean flow when the stratification is stable (respectively unstable) (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011).
Appendix B. Depth-independent variables
At a statistical steady state, stratified pressure-driven flows between solid boundaries have linearly varying shear stress 
$\bar {\tau }=\overline {\tau _{\nu }+\tau _w}$, with 
$\tau _{\nu }=\overline {\partial _z u}$ the viscous stress and 
$\tau _w=-\overline {wu}$ the Reynolds stress, and depth-independent heat flux 
$\bar {q}=\overline {wT-\partial _zT}$, where the overbar denotes horizontal and time averaging. These conservation equations for the vertical fluxes of momentum and heat are at the origin of the definitions of the friction velocity and Nusselt number, which typically read as 
$u_*=\sqrt {-\overline {(\tau _{\nu }+\tau _w)}}\Big|_{z=1}$ and 
$Nu=\int _0^{1} \bar {q} \,\textrm {d} z$ (assuming 
$z=1$ is the top of the fluid), respectively. With the phase-field method, these conservation equations are modified and some of the modifications are reflected in the definitions of 
$u_*$ and 
$Nu$ in (2.9a–c). In particular, 
$u_*$ in (2.9a–c) includes the linear damping term 
$\tau _d=-\int _0^{z} (1-\phi )u/\varGamma \,\textrm {d} z'/z$ that comes from the last term on the right-hand side of (2.5a). The true conservation of vertical momentum and heat fluxes based on governing equations (2.5) read as
\begin{equation} \partial_z \overline{\left(\tau + \tilde{\tau} \right)} = -2Pr^2Re, \quad \partial_z \overline{\left(q+\tilde{q} \right)} = 0 , \end{equation}
where 
$\tilde {\tau }$ and 
$\tilde {q}$ are the anomalous stress and heat flux due to damping of the advective terms in (2.5), i.e.
\begin{align} \tilde{\tau} = \int_0^{z} -(1-\phi)({\boldsymbol{u}}\cdot\boldsymbol{\nabla})u\,\textrm{d} z'/z,\quad \tilde{q} = \int_0^{z} -(1-\phi)({\boldsymbol{u}}\cdot\boldsymbol{\nabla})T\,\textrm{d} z'/z. \end{align}The existence of anomalous stress and heat fluxes means that the friction velocity and Nusselt numbers as defined in (2.9a–c) are based on a total stress and heat flux, which are not rigorously linearly varying or depth invariant.
We show in figure 15(a–c) the Reynolds stress 
$\overline {\tau _w}$, the viscous stress 
$\overline {\tau _{\nu }}$, the linear damping stress 
$\overline {\tau _d}$ and the Reynolds stress plus the anomalous stress 
$\overline {\tau _w+\tilde {\tau }}$ for stable, neutral and unstable stratification, respectively. Importantly, 
$\overline {\tau _w}$ and 
$\overline {\tau _w+\tilde {\tau }}$ overlap well, showing that the anomalous stress is negligible. The results of figure 15(d–f) further confirm that the anomalous stress is negligible in all simulations: the (approximate) total stress (solid lines) decreases linearly with 
$z$ in all stages and overlaps well with 
$\overline {\tau _w+\tau _{\nu }+\tau _d+\tilde {\tau }}$, i.e. the total stress (thin dashed lines) that includes the anomalous stress. We show in figure 15(g–i) 
$\bar {q}$ (solid lines) and 
$\overline {q+\tilde {q}}$ (thin dashed lines). For stable and neutral stratification (top and middle rows), 
$\bar {q}$ and 
$\overline {q+\tilde {q}}$ are constants with depth and overlap perfectly, suggesting that the anomalous heat flux is negligible. For unstable stratification (bottom row), we obtain similar results for stages I and II. For stage III, however, 
$\bar {q}$ is not perfectly constant, but deviates from 
$\overline {q+\tilde {q}}$ and peaks at 
$z\approx 1.15$, which is roughly the height of the channels. The relative discrepancy between 
$\bar {q}$ and 
$\overline {q+\tilde {q}}$ is of the order of 5 % and is a result of the damping of the advective terms in the momentum and heat equations (2.5). We expect that this discrepancy would decrease with increased resolution. Previous studies have alternatively considered advective terms with the same damped form as here, with a divergence damped form, i.e. 
$({\boldsymbol {u}}\cdot \boldsymbol {\nabla })(\phi {\boldsymbol {u}})$, or without any damping. There is no proof that any of these methods is more efficient than the other two. However, we would recommend using either one of the latter two methods, i.e. not the method used in this paper, in order to simplify the analysis of the shear stress and heat flux.
(a–c) Reynolds stress 
$\tau _w$ (dotted lines), viscous stress 
$\tau _{\nu }$ (dashed lines), linear damping stress 
$\tau _d$ (solid lines) and Reynolds stress plus anomalous stress 
$\tau _w+\tilde {\tau }$ (for stages II and III; black dashed lines) averaged in time over 20 friction time units and horizontal planes as functions of depth 
$z$ for stable, neutral and unstable stratification, respectively (see the text for more details). Blue, orange and green colours denote results obtained in stages I, II and III, respectively (same in (d–f) and (g–i)). (d–f) Total stress, i.e. 
$\tau _w+\tau _{\nu }+\tau _d$ (solid lines), and total stress plus the anomalous stress, i.e. 
$\tau _w+\tau _{\nu }+\tau _d+\tilde {\tau }$ (thin black dashed lines), averaged in time over 20 friction time units and horizontal planes as functions of depth 
$z$ for stable, neutral and unstable stratification, respectively. Note that the full stress is shifted to the right by 10 000 (20 000) between stage II and stage I and between stage III and stage II for the case of stable (unstable) stratification for clarity as all curves overlap otherwise. (g–i) Heat flux 
$q$ (solid lines) and heat flux plus anomalous heat flux 
$q+\tilde {q}$ (thin black dashed lines) averaged in time over 20 friction time units and horizontal planes as functions of depth 
$z$ for stable, neutral and unstable stratification, respectively. Note that the heat flux is shifted to the right by 2.5 (5) between stage II and stage I and between stage III and stage II for the case of stable (unstable) stratification for clarity.
Appendix C. Additional details on the simulation stages
In this section we give additional details on the simulation stages and sub-stages. For 
$t\leq t_*^{Ic}$ (including stage I), we solve (2.5a), (2.5b) and (2.5d) with 
$\phi \equiv 1$ and a no-slip isothermal top-boundary condition, i.e. 
${\boldsymbol {u}}=\boldsymbol {0}$ and 
$T=0$ at 
$z=1$. We use a straightforward half-channel flow configuration, i.e. without a solid domain, with 64 Chebyshev modes in the vertical direction. In stage II we add a solid layer of thickness 0.5 on top of the fluid domain and we use a compound Chebyshev basis stitched at 
$z=1.2$ with 256 (respectively 32) Chebyshev modes in the lower (respectively upper) region. The compound Chebyshev basis allows us to have a high vertical resolution near the interface's initial position. We solve (2.5a), (2.5b) and (2.5d) with 
$\phi$ prescribed, i.e. not varying in time (cf. the main text). In stage III we solve (2.5) with all variables freely evolving and we use the same spectral resolution as in stage II.
Our simulations until 
$t = t_*^{Ic}$ can be broken down into three sub-stages. In stage Ia, i.e. for 
$t_* < t_*^{Ia}$ (cf. light grey colour in figure 2), we run a low-resolution (128 Fourier modes in the 
$x$ and 
$y$ directions and 32 Chebyshev modes in the 
$z$ direction) spin-up simulation of an initially laminar flow superposed with three-dimensional velocity perturbations and no buoyancy effects (
$Ri_*=0$). In stage Ib, i.e. for 
$t_*^{Ia}\leq t_* < t_*^{Ib}$, we increase the resolution (256 Fourier modes in the 
$x$ and 
$y$ directions and 64 Chebyshev modes in the 
$z$ direction) but keep 
$Ri_*=0$ (cf. dark grey colour in figure 2). In stage I we turn on buoyancy effects (
$t_*^{Ib}\leq t_* < t_*^{Ic}$; cf. blue colour in figure 2), i.e. we use
\begin{equation} Ri_*(t) = Ri_*\tanh\left(t_*-t_*^{Ib} \right), \end{equation}
such that 
$Ri_*$ transitions smoothly (over the time scale of one friction time unit) from 0 at the end of stage Ib to the target value listed in table 1. Note that in the case of stabilizing buoyancy effects, we found that the turbulent flow relaxes to a laminar state when using (C1). In order to avoid this we used an intermediate stage with a more moderate target 
$Ri_*=20$ (sub-stage Ici) before transitioning to 
$Ri_*=40$, using a similar equation as (C1).
Appendix D. Higher-order interface evolution model
In this section we derive a reduced model for the evolution of the mean interface position 
$\bar {\xi }$, or ice thickness 
$h=h_0-(\bar {\xi }-{\xi }_0)$, with 
$h_0=1/2$ and 
${\xi }_0=1$ the initial ice thickness and interface position, which takes into account interface deformation. Under the steady-state assumption (instantaneous temperature diffusion), the evolution of 
$\bar {\xi }$ is controlled by the difference between the input heat flux from the fluid, 
$q^f$, and the mean heat flux at the ice top, 
$q_{top}$, i.e.
\begin{equation} \frac{\textrm{d}\bar{\xi}}{\textrm{d} t} = q^f-q_{top}. \end{equation}
At the leading order we assume that the interface is flat such that 
$q_{top}=h_0q^s/h$, which yields (3.4). At higher order, we take into account interface deformation, which can change 
$q_{top}$, using regular perturbation (Favier et al. Reference Favier, Purseed and Duchemin2019). Specifically, we seek a solution of the 
$x$-averaged steady-state heat equation, i.e. (dropping the tilde for 
$x$-averaged variables)
\begin{gather} \nabla^2 T = 0,\quad \xi(y) \leq z \leq L_z, \end{gather}
\begin{gather}T = -h_0q^s,\quad z=L_z, \end{gather}
\begin{gather}T = 0,\quad z=\xi(y), \end{gather}using a perturbation series of the form
\begin{equation} T(y,z) = T^{(0)}(z) + T^{(1)}(y,z) + T^{(2)}(y,z) + \cdots, \end{equation}
with 
$\partial _yT^{(0)}\equiv 0$ and 
$T^{(i)} \sim O(\epsilon ^i)$, where 
$\epsilon \ll 1$ is the dimensionless amplitude of the interface deformation. For simplicity, here we approximate the interface deformation as
\begin{equation} \xi(y) = \bar{\xi} + \epsilon \cos ky. \end{equation}The leading-order solution is
\begin{equation} T^{(0)} = -\left(\frac{z-\bar{\xi}}{L_z-\bar{\xi}}\right)h_0q^s, \end{equation}the first-order solution is
\begin{equation} T^{(1)} = \left(\frac{\epsilon h_0q^s}{L_z-\bar{\xi}} \right)\cos kx \frac{\sinh k\left(L_z-z \right)}{\sinh k\left(L_z-\bar{\xi}\right)}, \end{equation}and the second-order solution is
\begin{equation} T^{(2)} = \frac{\epsilon^2h_0q^s}{2\left(L_z-\bar{\xi}\right)\tanh k \left(L_z-\bar{\xi} \right)} \left[ \left(\frac{L_z-z}{L_z-\bar{\xi}}\right)+ \cos 2kx \frac{\sinh 2k\left(L_z-z \right)}{\sinh 2k\left(L_z-\bar{\xi}\right)} \right]. \end{equation}Thus, the second-order-accurate formula for the mean heat flux at the top of the ice reads as
\begin{equation} q_{top} = -\overline{\partial_z T}(z=L_z) \approx \frac{h_0q^s}{h}\left[ 1+ \frac{\epsilon^2k}{2h\tanh kh} \right], \end{equation}
which differs from the leading-order heat flux only at second order. The prediction for the evolution of the mean interface position for unstable stratification with 
$q^f=Nu^{III}$, which is shown by the red dash–dot lines in figure 11(a), is based on (D1) with 
$q_{top}$ given by (D10) and with 
$\epsilon =0.137$ and 
$k=2$ (as obtained from best-fit of the true interface topography at steady state for unstable stratification). We note that the quasi-steady-state assumption may affect the prediction of the transient evolution of the mean interface position adversely but has no effect on the final value, which is the primary goal of the reduced model.
Open access
