# Topography generation by melting and freezing in a turbulent shear flow

Louis-Alexandre Couston, Eric Hester, Benjamin Favier, John R Taylor, Paul R Holland, Adrian Jenkins

► To cite this version:

Louis-Alexandre Couston, Eric Hester, Benjamin Favier, John R Taylor, Paul R Holland, et al. Topography generation by melting and freezing in a turbulent shear flow. Journal of Fluid Mechanics, 2021, 911, 10.1017/jfm.2020.1064. hal-03375670

HAL Id: hal-03375670

https://hal.science/hal-03375670

Submitted on 13 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Topography generation by melting and freezing in a turbulent shear flow

Louis-Alexandre Couston<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,†, Eric Hester<fnref n="4" />, Benjamin Favier<fnref n="5" />, John R. Taylor<fnref n="2" />, Paul R. Holland<fnref n="1" /> and Adrian Jenkins<fnref n="1" />,<fnref n="6" />

(Received 20 April 2020; revised 23 November 2020; accepted 24 November 2020)

We report an idealized numerical study of a melting and freezing solid adjacent to a turbulent, buoyancy-affected shear flow, in order to improve our understanding of topography generation by phase changes in the environment. We use the phase-field method to dynamically couple the heat equation for the solid with the Navier–Stokes equations for the fluid. We investigate the evolution of an initially flat and horizontal solid boundary overlying a pressure-driven turbulent flow. We assume a linear equation of state for the fluid and change the sign of the thermal expansion coefficient, such that the background density stratification is either stable, neutral or unstable. We find that channels aligned with the direction of the mean flow are generated spontaneously by phase changes at the fluid–solid interface. Streamwise vortices in the fluid, the interface topography and the temperature field in the solid influence each other and adjust until a statistical steady state is obtained. The crest-to-trough amplitude of the channels is larger than approximately $10\delta_v$ in all cases, with $\delta_v$ the viscous length scale, but is much larger and more persistent for an unstable stratification than for a neutral or stable stratification. This happens because a stable stratification makes the cool melt fluid buoyant such that it shields the channel from further melting, whereas an unstable stratification makes the cool melt fluid sink, inducing further melting by rising hot plumes. The statistics of flow

† Email address for correspondence: louis.couston@ens-lyon.fr

<footnote n="1">¹British Antarctic Survey, Cambridge CB3 0ET, UK</footnote>

<footnote n="2">²Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, UK</footnote>

<footnote n="3">³Univ Lyon, Ens de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, F-69342 Lyon, France</footnote>

<footnote n="4">⁴Department of Mathematics and Statistics, University of Sydney, Sydney, New South Wales 2006, Australia</footnote>

<footnote n="5">⁵Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE, Marseille, France</footnote>

<footnote n="6">⁶Department of Geography and Environmental Sciences, Northumbria University, Newcastle upon Tyne NE1 8ST, UK</footnote>
velocities and melt rates are investigated, and we find that channels and keels emerging in
our simulations do not significantly change the mean drag coefficient.

**Key words:** geophysical and geological flows, turbulent flows, solidification/melting

**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. 2012; Rignot et al. 2013; Kennicutt 2019), 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 2019). 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. 2016). 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. 2019), but an ice shelf around Antarctica typically melts at a rate of only a few centimetres per day or less (Dutrieux et al. 2014). 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 2019), 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. 2013; Gourmelen et al. 2017) and terraces (Dutrieux et al.
2014) 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 2011), and, more generally,
rough features can be seen from the centimetre scale to tens-of-metre scale under sea
ice (Wadhams, Wilkinson & McPhail 2006; McPhee 2008; Lucieer et al. 2016) and up
to the kilometre scale under ice shelves (Nicholls 2006). 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. 2014) and recent laboratory experiments on
ice scallops (Bushuk et al. 2019).

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°C < T < 4°C at
atmospheric pressure (Thoma et al. 2010), 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).
Melting and freezing in a turbulent shear flow

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. 2018; Favier, Purseed & Duchemin 2019). 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 1994; Sullivan, Liu & Ecke 1996) and numerical simulations (Philippi et al. 2019). Streamwise patterns also emerge in dissolution experiments when the phase boundary is not perpendicular to gravity but inclined (Allen 1971; Cohen et al. 2020).

Despite the existence of many studies on pressure-driven and shear flows (Kelly 1994; Zonta & Soldati 2018), 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 (1980) 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 2019). However, Bushuk et al. (2019), who conducted an experiment similar to Gilpin et al. (1980), 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° to +90° compared with the maximum (crest) of boundary topography (Hanratty 1981). Recently, Claudin, Durán & Andreotti (2017) 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. (2017), 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.* 1993; Karma & Rappel 1998; Beckermann *et al.* 1999), but which was recently applied to the case of vigorous convective flows (Favier *et al.* 2019; Purseed *et al.* 2020). Other methods that simultaneously solve for the evolution of a fluid phase and a solid phase include the enthalpy method (Ulvrová *et al.* 2012), the level-set method (Gibou *et al.* 2007), the lattice-Boltzmann method (Rabbanipour Esfahani *et al.* 2018) and two-domain moving-boundary methods (Ulvrová *et al.* 2012). 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.* (2017) and experimentally by Bushuk *et al.* (2019). 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.* 2017; Bushuk *et al.* 2019).

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*<sub>x</sub> = 4π*H* and the transverse width is *L*<sub>y</sub> = 2π*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 (<sup>*f*</sup>) and (<sup>*s*</sup>) to denote variables in the fluid and the solid, respectively. The fluid velocity *u*(<sup>*f*</sup>) and pressure *p*(<sup>*f*</sup>) 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 *ρ*<sub>f</sub>, the same specific heat capacity at constant pressure *c*<sub>p</sub> and the same thermal conductivity *k*. Thus, the temperatures *T*(<sup>*f*</sup>) and *T*(<sup>*s*</sup>) evolve
<figure><img src="image_3.png" /><figcaption>Figure 1. (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\pi$, $l_y = 2\pi$ and $l_z = 1.5$, respectively.</figcaption></figure>

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

$$ \rho^{(f)} = \rho_f \left(1 - \alpha T^{(f)}\right), \qquad (2.1) $$

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.

$$ T^{(f)} = T^{(s)} = T_m, \qquad (2.2a) $$

$$ \rho_s \mathcal{L} v_n = q_n^{(f)} - q_n^{(s)} = -k_f \hat{n} \cdot \nabla T^{(f)} + k_s \hat{n} \cdot \nabla T^{(s)}, \qquad (2.2b) $$

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{n}$), $\mathcal{L}$ is the latent heat of fusion per unit mass, $q_n$ is the heat flux in direction $\hat{n}$, $k_s$ (respectively $k_f$) is the thermal conductivity of the solid (respectively fluid) and $\nabla$ is the gradient operator (Worster 2000). 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 \text{ kg m}^{-3}$, $c_{pf} \approx 4200 \text{ J kg}^{-1} \text{ K}^{-1}$ and $k_f \approx 0.6 \text{ W m}^{-1} \text{ K}^{-1}$, while $\rho_s \approx 917 \text{ kg m}^{-3}$, $c_{ps} \approx 2100 \text{ J kg}^{-1} \text{ K}^{-1}$ and $k_s \approx 2.2 \text{ W m}^{-1} \text{ 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.

Here we use a volume-penalization method (Angot, Bruneau & Fabrie 1999), 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.

$$ \frac{\partial \mathbf{u}}{\partial t} + \phi (\mathbf{u} \cdot \nabla) \mathbf{u} = \nu \nabla^2 \mathbf{u} - \frac{1}{\rho_f} \nabla p + \alpha g T \hat{\mathbf{z}} + \frac{\Pi}{\rho_f} \hat{\mathbf{x}} - \frac{(1-\phi)}{\tau_p} \mathbf{u}, \quad (2.3a) $$

$$ \frac{\partial T}{\partial t} + \phi (\mathbf{u} \cdot \nabla) T = \kappa \nabla^2 T - \frac{\mathcal{L}}{c_p} \frac{\partial \phi}{\partial t}, \quad (2.3b) $$

$$ \frac{\partial \phi}{\partial t} = a \nabla^2 \phi + b \phi (1 - \phi) [2\phi - 1 + c(T - T_m)], \quad (2.3c) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (2.3d) $$

with $\mathbf{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), $v$ is the constant kinematic viscosity, $g$ is the standard gravity, $\Pi$ 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{\mathbf{z}}$ and $\hat{\mathbf{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. 2019), 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

$$ \left. \begin{array}{l} (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 + \Delta T \tilde{T}, \quad p = \rho_f u_\kappa^2 \tilde{p}, \quad \phi = \tilde{\phi}, \end{array} \right\} (2.4a-f) $$

with $u_\kappa = \kappa/H$ the diffusion velocity scale, and $\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

$$ \frac{\partial \mathbf{u}}{\partial t} + \phi (\mathbf{u} \cdot \nabla) \mathbf{u} = Pr \nabla^2 \mathbf{u} - \nabla p + PrRaT\hat{z} + 2Pr^2Re\hat{x} - Pr \frac{(1-\phi)}{\Gamma} \mathbf{u}, \quad (2.5a) $$

$$ \frac{\partial T}{\partial t} + \phi (\mathbf{u} \cdot \nabla) T = \nabla^2 T - St \frac{\partial \phi}{\partial t}, \quad (2.5b) $$

$$ \frac{\partial \phi}{\partial t} = A \nabla^2 \phi + B\phi(1-\phi)(2\phi - 1 + CT), \quad (2.5c) $$

$$ \nabla \cdot \mathbf{u} = 0. \quad (2.5d) $$

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

$$ Pr = \frac{\nu}{\kappa}, \quad Re = \frac{\Pi H^3}{2\rho_f \nu^2}, \quad Ra = \frac{\alpha g \Delta T H^3}{\nu \kappa}, \quad St = \frac{\mathcal{L}}{c_p \Delta T}. \quad (2.6a-d) $$

The additional parameters $\Gamma = \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. $\mathbf{u} = \mathbf{0}$ and $T = T_t < 0$ at $z = 1.5$. We impose free-slip, fixed-temperature conditions on the bottom boundary, i.e. $\partial_z \mathbf{u} = \partial_z \nu = 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\pi, 2\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.

$$ Re_* = \sqrt{2Re} = \sqrt{\frac{\Pi H^3}{\rho_f \nu^2}}, \quad Ri_* = \frac{-Ra}{PrRe} = \frac{-2\rho_f \alpha \Delta T_g}{\Pi}, \quad (2.7a,d) $$

since they are more commonly used than Re and Ra in turbulent channel flow studies (García-Villalba & del Álamo 2011; Zonta & Soldati 2018). 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 $[\Pi 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{\Pi 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

$$L_{MO} = \frac{Re_*^3 Pr^2}{Nu Ra}, \qquad (2.8)$$

in terms of dimensionless variables and which is often reported in mixed-convection experiments (Pirozzoli et al. 2017; Blass et al. 2020), 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 = 11250$), $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 1961; Gage & Reid 1968). 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. 2020). 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 $\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 (1998) 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-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
<table><caption>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$.</caption><thead><tr><th>Stratification</th><th>Ri<sub>*</sub></th><th>Re<sub>b</sub><sup>I</sup></th><th>Re<sub>b</sub><sup>III</sup></th><th>Nu<sup>I</sup></th><th>q<sup>s</sup></th><th>Nu<sup>III</sup></th><th>ξ<sup>III</sup></th><th>(ξ<sub>−</sub>/δ<sub>v</sub>, ξ<sub>+</sub>/δ<sub>v</sub>)</th><th>10<sup>3</sup>C<sub>D</sub><sup>I</sup></th><th>10<sup>3</sup>C<sub>D</sub><sup>III</sup></th></tr></thead><tbody><tr><td>Stable</td><td>40</td><td>2720</td><td>2630</td><td>3.02</td><td>2.83</td><td>3.09</td><td>1.042</td><td>(2.4,2.2)</td><td>5.9</td><td>6.7</td></tr><tr><td>Neutral</td><td>0</td><td>2290</td><td>2240</td><td>4.58</td><td>4.48</td><td>4.75</td><td>1.025</td><td>(5.2,4.1)</td><td>8.6</td><td>9.2</td></tr><tr><td>Unstable</td><td>−40</td><td>1970</td><td>1910</td><td>6.74</td><td>6.69</td><td>7.75</td><td>1.037</td><td>(30,15)</td><td>11.3</td><td>12.7</td></tr></tbody></table>

over a very thin diffuse interface of thickness $\approx 0.007$. Simulation parameters and output variables are provided in table 1.

## 2.3. Variables of interest

We define the friction velocity, the bulk velocity and the Nusselt number as

$$u_* = \langle \phi \rangle \sqrt{-\overline{(\tau_d + \tau_v + \tau_w)}} \Big|_{z=1.5}, \quad u_b = \frac{\langle u \rangle}{\langle \phi \rangle}, \quad Nu = \langle q \rangle, \quad (2.9a-c)$$

respectively (cf. details in appendix B), where the overbar denotes the horizontal average and $\langle \cdot \rangle = \int_V dV/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_v$ and $\tau_w$ are the linear damping, viscous and Reynolds shear stresses, and $q = wT - \partial_z T$ is the heat flux. At statistical steady state, the full shear stress $\bar{\tau} = (\tau_d + \tau_v + \tau_w)$ 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 \le 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

$$\xi(x, y) = \int_{0}^{l_z} \phi \, dz, \quad (2.10)$$

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.

$$C_D = 2 \left( \frac{u_*}{u_b} \right)^2 \quad (2.11)$$

(García-Villalba & del Álamo 2011; Pirozzoli et al. 2017). 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_*^{lb} \le t \le t_*^{lc}$) 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_*^{lc} < t \le 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 \le z \le 1.5$ and the phase field is prescribed as $\phi = 0.5 \{1 - \tanh[2(z-1)/\delta]\}$, 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

$$T = -q^s(z-1), \quad 1 \le z \le 1.5, \qquad (3.1)$$

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

$$ \left. \begin{aligned} Re_b^I &= \int_{t_*^c - \Delta_*}^{t_*^c} \frac{u_b \, dt_*}{Pr \Delta_*}, & Re_b^{III} &= \int_{t_*^{III} - \Delta_*}^{t_*^{III}} \frac{u_b \, dt_*}{Pr \Delta_*}, \\ Nu^I &= \int_{t_*^c - \Delta_*}^{t_*^c} \frac{Nu \, dt_*}{\Delta_*}, & Nu^{III} &= \int_{t_*^{III} - \Delta_*}^{t_*^{III}} \frac{Nu \, dt_*}{\Delta_*} \end{aligned} \right\} (3.2a-d) $$
Melting and freezing in a turbulent shear flow

<figure><img src="image_4.png" /><figcaption>Figure 2. 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_*^b, t_*^{ci}]$ in (a) during which $Ri_* = 20$. See the text and appendix C for more details.</figcaption></figure>
with $\Delta_* = 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.

Buoyancy effects are turned off for $t_* \le t_*^{lb}$, such that the results of figure 2 are exactly the same for all three simulations until $t_* = t_*^{lb}$. Upon turning on buoyancy, i.e. for $t_* \ge t_*^{lb}$ (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 2011; Pirozzoli et al. 2017), and the important point is that the heat flux is the variable that changes the most with buoyancy effects. Here, $Nu^l = 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_*^{lc}$ 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 1999). The drag coefficient is of order $10^{-2}$ and decreases (respectively increases) significantly at $t_* = t_*^{lb}$, 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 2011; Pirozzoli et al. 2017). 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).
Melting and freezing in a turbulent shear flow

<figure><img src="image_5.png" /><figcaption>Figure 3. 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.</figcaption></figure>

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 ($\tilde{\cdot}$) and variables averaged in the $x$ direction minus the horizontal mean by a prime ($'$), such that e.g. $\xi' = \tilde{\xi} - \tilde{\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_v = 1/Re_*$ for stable, neutral and unstable stratification, respectively (note that $\delta_v$ 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_v$, 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).

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
<figure><img src="image_6.png" /><figcaption>Figure 4. (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} - \tilde{\xi}_{\infty}$ normalized by the viscous length scale $\delta_v = 1/Re_*$ and as functions of $t_*$.
The top, middle and bottom rows show the results for stable, neutral and unstable stratification, respectively.</figcaption></figure>

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
# Melting and freezing in a turbulent shear flow

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} = \tilde{w}\tilde{T} - \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{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,
<figure><img src="image_7.png" /><figcaption>Figure 5. 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_v = 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.</figcaption></figure>
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.

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. 2017; Zonta & Soldati 2018), but their amplitude clearly varies with $R_i$. 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. 2017), 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_* \ge 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.

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_y k_y / (2\pi) = 2$,
<figure><img src="image_8.png" /><figcaption>Figure 6. (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{m}$ averaged in the streamwise direction as a function of $(t_*, y)$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. 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.</figcaption></figure>

which is close to the critical wavenumber $L_y k_c / (2\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
Melting and freezing in a turbulent shear flow

<figure><img src="image_10.png" /><figcaption>Figure 8. 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.</figcaption></figure>

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.
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. For 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.

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
<figure><img src="image_11.png" /><figcaption>Figure 9. Probability density functions of the streamwise velocity fluctuations $u - \bar{u}$ (leftmost panels), vertical velocity fluctuations $w - \bar{w}$ (middle panels) and temperature gradient fluctuations $\partial_z T - \bar{\partial_z T}$ (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{u}$ and $w - \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_z T - \bar{\partial_z T}$ 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{u}$ and $w - \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.</figcaption></figure>

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

$$Ri_*(t_*) = -40[1 - f(t_*, 5)] + 20[2 + f(t_*, 12) + f(t_*, 19)], \quad (3.3)$$

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. 2017). 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 2011), while the decreasing topography amplitude is expected to result in flattening and weakening vortices. Ultimately, the topography disappears and the vortices weaken significantly.

### 3.5. Melting

In this section we discuss the evolution of the mean interface position $\tilde{\xi}$ with time and the statistics of melting $\dot{m} = \partial_t \tilde{\xi}$ at the statistical steady state. We first show the evolution
Melting and freezing in a turbulent shear flow

<figure><img src="image_12.png" /><figcaption>Figure 10. 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}$.</figcaption></figure>

of $\tilde{\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 $\tilde{\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^l$) 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^l - 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 $\tilde{\xi}$ is faster for stable than for unstable stratification.
<figure><img src="image_13.png" /><figcaption>Figure 11. (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 $\Delta_* = 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$.</figcaption></figure>

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.

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

$$ \frac{d\bar{\xi}}{dt} = q^f - \frac{h_0 q^s}{h}, \qquad (3.4) $$

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^*$. 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^* + 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)
Melting and freezing in a turbulent shear flow

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 $\tilde{\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 $\tilde{\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.

The p.d.f.s for the melt rate $\dot{m}$ are shown in figure 12(b). The temporal and spatial average, $\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
<figure><img src="image_14.png" /><figcaption>Figure 12. 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).</figcaption></figure>

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. 2010). 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 2016).

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°C. If we further assume that $H = 10$ cm and $v = \kappa = 10^{-6} \text{ m}^2 \text{ 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 \text{ cm s}^{-1}$ and the thermal expansion coefficient is $5 \times 10^{-7} \text{ 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. 2017), 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. 2017; Blass et al. 2020), 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 \ge 0.1$, and to grow in size as $Ri_b$ increases. In the limit $Ri_b \to \infty$, channels and keels may be expected to eventually disappear. Indeed, as $Ri_b \to \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. 2018). 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 2015), 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_* \ge 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_* \ge 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. 2019). 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 \ge O(1000-10000)$ and have always reported a scallop wavelength smaller than the fluid depth, i.e. $\lambda < H$ (Blumberg & Curl 1974; Thomas 1979; Claudin et al. 2017). Considering the upper limit $\lambda = H$ means that scallops are predicted to emerge for $Re_* \ge O(1000-10000)$ 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. 2019), in coupled fluid-solid simulations would be a major achievement, which could complement fluid-only simulations at planar ice boundaries (Gayen, Griffiths & Kerr 2016; Keitzl, Mellado & Notz 2016a,b; Mondal et al. 2019; Vreugdenhil & Taylor 2019). 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³ 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
Melting and freezing in a turbulent shear flow

verify or refine the most recent theoretical predictions on scallop formation and saturation (Claudin *et al.* 2017) 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*<sub>*</sub> 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.* 2017). 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 (1971)
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 1994; Pabiou, Mergui & Bénard 2005; Blass et al. 2020). 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 2013; Courrech du Pont, Narteau & Gao 2014).

**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 Skłodowska-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.

Author ORCIDs.

D Louis-Alexandre Couston https://orcid.org/0000-0002-2184-2472;

D Benjamin Favier https://orcid.org/0000-0002-1184-2989;

D John R. Taylor https://orcid.org/0000-0002-1292-3756;

D Adrian Jenkins https://orcid.org/0000-0002-9117-0616.

**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 *Γ* of the phase-field equations (2.5) must satisfy several constraints (Hester *et al.* 2020), 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
L.-A. Couston and others

<figure><img src="image_15.png" /><figcaption>Figure 13. 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.</figcaption></figure>

$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, $\Gamma = (\delta/2.648228)^2$ and we require time steps to be always smaller than $\Gamma/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.

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_v = 1/Re_*$, Here, we chose to have $\delta \approx \delta_v$ in all simulations, i.e. $\delta$ is equal to 1 wall unit $\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_v/2$. The damping time scale $\Gamma$ in (2.5) is set to $\Gamma = (\delta/2.648228)^2$ in order to cancel first-order errors in the phase-field model (Hester et al. 2020).

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^+ \ge 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
<figure><img src="image_16.png" /><figcaption>Figure 14. Wall-normalized velocity $U^+ = Pru/u_*$ (circles, left axes) and temperature $T^+ = u_*(T/(PrNu))$ (crosses, right axes) expressed in terms of wall units as functions of the wall coordinate $z^+ = 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_*^b, t_*^c]$) whereas orange symbols denote results obtained with volume penalization (i.e. for $t_* \in [t_*^c, t_*^h]$). 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.</figcaption></figure>

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 2011).

## Appendix B. Depth-independent variables

At a statistical steady state, stratified pressure-driven flows between solid boundaries have linearly varying shear stress $\bar{\tau} = \bar{\tau}_v + \bar{\tau}_w$, with $\bar{\tau}_v = \partial_z u$ the viscous stress and $\bar{\tau}_w = -\overline{w\bar{u}}$ the Reynolds stress, and depth-independent heat flux $\bar{q} = \overline{wT - \partial_z T}$, 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{-(\bar{\tau}_v + \bar{\tau}_w)}|_{z=1}$ and $Nu = \int_0^1 \bar{q} dz$ (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 $\bar{\tau}_d = -\int_0^z (1-\phi)u/\Gamma dz'/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

$$ \partial_z(\bar{\tau} + \bar{\tau}) = -2Pr^2Re, \quad \partial_z(q + \bar{q}) = 0, \qquad (B1a,b) $$
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.

$$ \tilde{\tau} = \int_{0}^{z} -(1-\phi)(\mathbf{u} \cdot \nabla)\mathbf{u} \, dz'/z, \quad \tilde{q} = \int_{0}^{z} -(1-\phi)(\mathbf{u} \cdot \nabla)T \, dz'/z. \quad (B2a,b) $$

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}_v$, 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 + \overline{\tau}_v + \overline{\tau}_d + \tilde{\tau}$, i.e. the total stress (thin dashed lines) that includes the anomalous stress. We show in figure 15(g-i) $\tilde{q}$ (solid lines) and $\overline{q} + \tilde{q}$ (thin dashed lines). For stable and neutral stratification (top and middle rows), $\tilde{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, $\tilde{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 $\tilde{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. $(\mathbf{u} \cdot \nabla)(\phi \mathbf{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.

## Appendix C. Additional details on the simulation stages

In this section we give additional details on the simulation stages and sub-stages. For $t \le t^c$ (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. $\mathbf{u} = \mathbf{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^c$ can be broken down into three sub-stages. In stage Ia, i.e. for $t_* < t^a$ (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_* \le t^a < t^b$, we increase the resolution
Melting and freezing in a turbulent shear flow

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Figure 15. (a–c) Reynolds stress $\tau_w$ (dotted lines), viscous stress $\tau_v$ (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_v + \tau_d$ (solid lines), and total stress plus the anomalous stress, i.e. $\tau_w + \tau_v + \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.</figcaption></figure>
(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_*^{lb} \le t_* < t_*^{lc}$; cf. blue colour in figure 2), i.e. we use

$$Ri_*(t) = Ri_* \tanh(t_* - t_*^{lb}), \qquad (C1)$$

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.

$$\frac{d\bar{\xi}}{dt} = q^f - q_{top}. \qquad (D1)$$

At the leading order we assume that the interface is flat such that $q_{top} = h_0q^f/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. 2019). Specifically, we seek a solution of the x-averaged steady-state heat equation, i.e. (dropping the tilde for x-averaged variables)

$$\nabla^2 T = 0, \quad \xi(y) \le z \le L_z, \qquad (D2a,b)$$

$$T = -h_0 q^s, \quad z = L_z, \qquad (D3a,b)$$

$$T = 0, \quad z = \xi(y), \qquad (D4a,b)$$

using a perturbation series of the form

$$T(y, z) = T^{(0)}(z) + T^{(1)}(y, z) + T^{(2)}(y, z) + \dots, \qquad (D5)$$

with $\partial_y T^{(0)} = 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

$$\xi(y) = \bar{\xi} + \epsilon \cos ky. \qquad (D6)$$

The leading-order solution is

$$T^{(0)} = - \left( \frac{z - \bar{\xi}}{L_z - \bar{\xi}} \right) h_0 q^s, \qquad (D7)$$

the first-order solution is

$$T^{(1)} = \left( \frac{\epsilon h_0 q^s}{L_z - \bar{\xi}} \right) \cos kx \frac{\sinh k(L_z - z)}{\sinh k(L_z - \bar{\xi})}, \qquad (D8)$$
and the second-order solution is

$$ T^{(2)} = \frac{\epsilon^2 h_0 q^s}{2 (L_z - \bar{\xi}) \tanh k (L_z - \bar{\xi})} \left[ \left( \frac{L_z - z}{L_z - \bar{\xi}} \right) + \cos 2kx \frac{\sinh 2k (L_z - z)}{\sinh 2k (L_z - \bar{\xi})} \right]. \quad (D9) $$

Thus, the second-order-accurate formula for the mean heat flux at the top of the ice reads as

$$ q_{top} = -\bar{\partial}_z T(z = L_z) \approx \frac{h_0 q^s}{h} \left[ 1 + \frac{\epsilon^2 k}{2h \tanh kh} \right], \quad (D10) $$

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.

REFERENCES

ALLEN, J.R.L. 1971 Bed forms due to mass transfer in turbulent flows: a kaleidoscope of phenomena. *J. Fluid Mech.* **49** (1), 49–63.

ANGOT, P., BRUNEAU, C.H. & FABRIE, P. 1999 A penalization method to take into account obstacles in incompressible viscous flows. *Numer. Math.* **81** (4), 497–520.

BECKERMANN, C., DIEPERS, H.J., STEINBACH, I., KARMA, A. & TONG, X. 1999 Modeling melt convection in phase-field simulations of solidification. *J. Comput. Phys.* **154** (2), 468–496.

BLASS, A., ZHU, X., VERZICCO, R., LOHSE, D. & STEVENS, R.J.A.M. 2020 Flow organization and heat transfer in turbulent wall sheared thermal convection. *J. Fluid Mech.* **897**, A22.

BLUMBERG, P.N. & CURL, R.L. 1974 Experimental and theoretical studies of dissolution roughness. *J. Fluid Mech.* **65** (4), 735–751.

BURNS, K.J., VASIL, G.M., OISHI, J.S., LECOANET, D. & BROWN, B.P. 2020 Dedalus: a flexible framework for numerical simulations with spectral methods. *Phys. Rev. Res.* **2** (2), 23068.

BUSHUK, M., HOLLAND, D.M., STANTON, T.P., STERN, A. & GRAY, C. 2019 Ice scallops: a laboratory investigation of the ice-water interface. *J. Fluid Mech.* **873**, 942–976.

CHANDRASEKHAR, S. 1961 *Hydrodynamic and Hydromagnetic Stability*. Dover.

CHARRU, F., ANDREOTTI, B. & CLAUDIN, P. 2013 Sand ripples and dunes. *Annu. Rev. Fluid Mech.* **45** (1), 469–493.

CLAUDIN, P., DURÁN, O. & ANDREOTTI, B. 2017 Dissolution instability and roughening transition. *J. Fluid Mech.* **832**, 832R21–832R214.

COHEN, C., BERHANU, M., DERR, J. & DU PONT, S. 2020 Buoyancy-driven dissolution of inclined blocks: erosion rate and pattern formation. *Phys. Rev. Fluids* **5** (5), 053802.

COURRECH DU PONT, S., NARTEAU, C. & GAO, X. 2014 Two modes for dune orientation. *Geology* **42** (9), 743–746.

DALLASTON, M.C., HEWITT, I.J. & WELLS, A.J. 2015 Channelization of plumes beneath ice shelves. *J. Fluid Mech.* **785**, 109–134.

DAVIS, P.E.D. & NICHOLLS, K.W. 2019 Turbulence observations beneath Larsen C ice shelf, Antarctica. *J. Geophys. Res.* **124** (8), 5529–5550.

DINNIMAN, M.S., ASAY-DAVIS, X.S., GALTON-FENZI, B.K., HOLLAND, P.R., JENKINS, A. & TIMMERMANN, R. 2016 Modeling ice shelf/ocean interaction in Antarctica: a review. *Oceanography* **29** (4), 144–153.

DUTRIEUX, P., STEWART, C., JENKINS, A., NICHOLLS, K.W., CORR, H.F.J., RIGNOT, E. & STEFFEN, K. 2014 Basal terraces on melting ice shelves. *Geophys. Res. Lett.* **41** (15), 5506–5513.

FAVIER, B., PURSEED, J. & DUCHEMIN, L. 2019 Rayleigh-Bénard convection with a melting boundary. *J. Fluid Mech.* **858**, 437–473.

GAGE, K.S. & REID, W.H. 1968 The stability of thermally stratified plane Poiseuille flow. *J. Fluid Mech.* **33** (1), 21–32.
GARCÍA-VILLALBA, M. & DEL ÁLAMO, J.C. 2011 Turbulence modification by stable stratification in channel flow. *Phys. Fluids* **23** (4), 045104.

GAYEN, B., GRIFFITHS, R.W. & KERR, R.C. 2016 Simulation of convection at a vertical ice face dissolving into saline water. *J. Fluid Mech.* **798**, 284–298.

GIBOU, F., CHEN, L., NGUYEN, D. & BANERJEE, S. 2007 A level set based sharp interface method for the multiphase incompressible Navier–Stokes equations with phase change. *J. Comput. Phys.* **222** (2), 536–555.

GILPIN, R.R., HIRATA, T. & CHENG, K.C. 1980 Wave formation and heat transfer at an ice-water interface in the presence of a turbulent flow. *J. Fluid Mech.* **99** (3), 619–640.

GOURMELEN, N., et al. 2017 Channelized melting drives thinning under a rapidly melting Antarctic ice shelf. *Geophys. Res. Lett.* **44** (19), 9796–9804.

HANRATTY, T.J. 1981 Stability of surfaces that are dissolving or being formed by convective diffusion. *Annu. Rev. Fluid Mech.* **13** (1), 231–252.

HESTER, E.W., COUSTON, L.-A., FAVIER, B., BURNS, K.J. & VASIL, G.M. 2020 Improved phase-field models of melting and dissolution in multi-component flows. *Proc. R. Soc. Lond. A* **476** (2242), 20200508.

HOBSON, B.W., SHERMAN, A.D. & MCGILL, P.R. 2011 Imaging and sampling beneath free-drifting icebergs with a remotely operated vehicle. *Deep-Sea Res. II* **58** (11–12), 1311–1317.

HOLLAND, D.M. & JENKINS, A. 1999 Modeling thermodynamic ice-ocean interactions at the base of an ice shelf. *J. Phys. Oceanogr.* **29** (8 PART 1), 1787–1800.

JENKINS, A. 2016 A simple model of the ice shelf-ocean boundary layer and current. *J. Phys. Oceanogr.* **46** (6), 1785–1803.

KARMA, A. & RAPPEL, W.-J. 1998 Quantitative phase-field modeling of dendritic growth in two and three dimensions. *Phys. Rev. E* **57** (4), 4323–4349.

KEITZL, T., MELLADO, J.P. & NOTZ, D. 2016a Impact of thermally driven turbulence on the bottom melting of ice. *J. Phys. Oceanogr.* **46** (4), 1171–1187.

KEITZL, T., MELLADO, J.P. & NOTZ, D. 2016b Reconciling estimates of the ratio of heat and salt fluxes at the ice-ocean interface. *J. Geophys. Res.* **121**, 8419–8433.

KELLY, R.E. 1994 The onset and development of thermal convection in fully developed shear flows. *Adv. Appl. Mech.* **31**, 35–112.

KENNICUTT, M.C., et al. 2019 Sustained Antarctic research: a 21st century imperative. *One Earth* **1** (1), 95–113.

KERR, R.C. 1994 Dissolving driven by vigorous compositional convection. *J. Fluid Mech.* **280**, 287–302.

LUCIEER, V., NAU, A.W., FORREST, A.L. & HAWES, I. 2016 Fine-scale sea ice structure characterized using underwater acoustic methods. *Remote Sens.* **8** (10), 1–17.

MCPHEE, M. 2008 *Air-Ice-Ocean Interaction: Turbulent Ocean Boundary Layer Exchange Processes*. Springer.

MOIN, P. & MAHESH, K. 1998 Direct numerical simulation: a tool in turbulence research. *Annu. Rev. Fluid Mech.* **30** (1), 539–578.

MONDAL, M., GAYEN, B., GRIFFITHS, R.W. & KERR, R.C. 2019 Ablation of sloping ice faces into polar seawater. *J. Fluid Mech.* **863**, 545–571.

NICHOLLS, K.W., et al. 2006 Measurements beneath an Antarctic ice shelf using an autonomous underwater vehicle. *Geophys. Res. Lett.* **33** (8), 2–5.

PABIOU, H., MERGUI, S. & BÉNARD, C. 2005 Wavy secondary instability of longitudinal rolls in Rayleigh–Bénard-Poiseuille flows. *J. Fluid Mech.* **542**, 175–194.

PHILIPPI, J., BERHANU, M., DERR, J. & DU PONT, S. 2019 Solutal convection induced by dissolution. *Phys. Rev. Fluids* **4** (10), 103801.

PIROZZOLI, S., BERNARDINI, M., VERZICCO, R. & ORLANDI, P. 2017 Mixed convection in turbulent channels with unstable stratification. *J. Fluid Mech.* **821**, 482–516.

PRITCHARD, H.D., LIGTENBERG, S.R.M., FRICKER, H.A., VAUGHAN, D.G., VAN DEN BROEKE, M.R. & PADMAN, L. 2012 Antarctic ice-sheet loss driven by basal melting of ice shelves. *Nature* **484** (7395), 502–505.

PURSEED, J., FAVIER, B., DUCHEMIN, L. & HESTER, E.W. 2020 Bistability in Rayleigh–Bénard convection with a melting boundary. *Phys. Rev. Fluids* **5** (2), 023501.

RABBANIPOUR ESFAHANI, B., HIRATA, S.C., BERTI, S. & CALZAVARINI, E. 2018 Basal melting driven by turbulent thermal convection. *Phys. Rev. Fluids* **3** (5), 053501.

RIGNOT, E., JACOBS, S., MOuginot, J. & SCHEUCHL, B. 2013 Ice-shelf melting around Antarctica. *Science* **341** (6143), 266–270.
Melting and freezing in a turbulent shear flow

STANTON, T.P., SHAW, W.J., TRUFFER, M., CORR, H.F.J., PETERS, L.E., RIVERMAN, K.L., BINDSCHADLER, R., HOLLAND, D.M. & ANANDAKRISHNAN, S. 2013 Channelized ice melting in the ocean boundary layer beneath pine island glacier, Antarctica. *Science* **341** (6151), 1236–1239.

SULLIVAN, T.S., LIU, Y. & ECKE, R.E. 1996 Turbulent solutal convection and surface patterning in solid dissolution. *Phys. Rev. E* **54** (1), 486–495.

SUTHERLAND, D.A., JACKSON, R.H., KIENHOLZ, C., AMUNDSON, J.M., DRYER, W.P., DUNCAN, D., EIDAM, E.F., MOTYKA, R.J. & NASH, J.D. 2019 Direct observations of submarine melt and subsurface geometry at a tidewater glacier. *Science* **365** (6451), 369–374.

THOMA, M., GROSFELD, K., SMITH, A.M. & MAYER, C. 2010 A comment on the equation of state and the freezing point equation with respect to subglacial lake modelling. *Earth Planet. Sci. Lett.* **294** (1–2), 80–84.

THOMAS, R.M. 1979 Size of scallops and ripples formed by flowing water. *Nature* **277** (5694), 281–283.

TOPPALADODDI, S. & WETTLAUFER, J.S. 2019 The combined effects of shear and buoyancy on phase boundary stability. *J. Fluid Mech.* **868**, 648–665.

ULVROVÁ, M., LABROSSE, S., COLTICE, N., RÁBACK, P. & TACKLEY, P.J. 2012 Numerical modelling of convection interacting with a melting and solidification front: application to the thermal evolution of the basal magma ocean. *Phys. Earth Planet. Inter.* **206–207**, 51–66.

VREUGDENHIL, C.A. & TAYLOR, J.R. 2019 Stratification effects in the turbulent boundary layer beneath a melting ice shelf: insights from resolved large-eddy simulations. *J. Phys. Oceanogr.* **49** (7), 1905–1925.

WADHAMS, P., WILKINSON, J.P. & MCPHAIL, S.D. 2006 A new view of the underside of Arctic sea ice. *Geophys. Res. Lett.* **33** (4), 1–5.

WANG, S.-L., SEKERKA, R.F., WHEELER, A.A., MURRAY, B.T., CORIELL, S.R., BRAUN, R.J. & MCFADDEN, G.B. 1993 Thermodynamically-consistent phase-field models for solidification. *Physica D* **69** (1), 189–200.

WELLS, A.J., HITCHEN, J.R. & PARKINSON, J.R.G. 2019 Mushy-layer growth and convection, with application to sea ice. *Phil. Trans. R. Soc. Lond. A* **377** (2146), 2018.0165.

WORSTER, M.G. 2000 Solidification of fluids. In *Perspectives in Fluid Dynamics* (ed. G. Worster, K. Moffatt & G. Batchelor), pp. 393–446. Cambridge University Press.

ZONTA, F. & SOLDATI, A. 2018 Stably stratified wall-bounded turbulence. *Appl. Mech. Rev.* **70** (4), 040801.