Topography generation by melting and freezing in a turbulent shear flow

Abstract 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 _{\nu }$ in all cases, with $\delta _{\nu }$ 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 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.


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

Phase-field method
We investigate the generation of topography due to uneven melting and freezing at a fluid-solid interface. The solid is fixed and located above the fluid where a Poiseuille/channel flow develops due to an external pressure gradient (see figure 1). The initial thickness of the fluid (respectively solid) layer is H (respectively H/2), such that the channel full depth is 3H/2. The domain length (in the direction of the flow) is L x = 4πH and the transverse width is L y = 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 ( f ) and (s) to denote variables in the fluid and the solid, respectively. The fluid velocity u ( f ) and pressure p ( f ) evolve according to the Navier-Stokes equations under the Boussinesq approximation. For simplicity, we assume that the solid and fluid phases have the same thermodynamical and transport properties, i.e. the same reference density ρ f , the same specific heat capacity at constant pressure c p and the same thermal conductivity k. Thus, the temperatures T ( f ) and T (s) evolve according to the same advection-diffusion (energy) equation, which turns into the heat equation where the velocity is zero. We consider a generic linear equation of state for the fluid, i.e. not specific to water, with the density related to temperature through with α 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.
where ρ 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 vectorn), L is the latent heat of fusion per unit mass, q n is the heat flux in directionn, k s (respectively k f ) is the thermal conductivity of the solid (respectively fluid) and ∇ 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 ρ s = ρ 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. ρ f ≈ 999 kg m −3 , c pf ≈ 4200 J kg −1 K −1 and k f ≈ 0.6 W m −1 K −1 , while ρ s ≈ 917 kg m −3 , c ps ≈ 2100 J kg −1 K −1 and k s ≈ 2.2 W m −1 K −1 . The relative differences are small however, i.e. within a factor of four or less, such that we do not expect fundamental differences between the physics in our model and real processes involving water and ice in nature.
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 φ, i.e. (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), ν is the constant kinematic viscosity, g is the standard gravity, Π is the imposed pressure-gradient force and κ = k/ρ f /c p is the constant thermal diffusivity; τ p , a, b and c are parameters related to volume penalization and the phase-field method, which we define later, andẑ andx 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, φ, also known as the phase-field variable or order parameter, satisfies a forced diffusion equation (2.3c) with parameters tuned such that φ 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 φ 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 φ 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.

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 τ κ = H 2 /κ as a normalizing time scale, i.e. the dimensionless variables, denoted by tildes, are defined as (x, y, z) = (Hx, Hỹ, Hz), t = τ κt , u = u κũ , (2.4a-f ) with u κ = κ/H the diffusion velocity scale, and ΔT = T b − T m is the temperature scale with T b the dimensional temperature on the bottom boundary. The time scale τ κ 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 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 The additional parameters Γ = τ p νH 2 /κ 2 , A = a/κ, B = b/(κ/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. u = 0 and T = T t < 0 at z = 1.5. We impose free-slip, fixed-temperature conditions on the bottom boundary, i.e. ∂ z u = ∂ z v = 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π, 2π, 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.
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 [Π H 2 (1 − z 2 /H 2 )/(2ρ f ν)]| z=0 using dimensional variables, while the latter is based on the friction velocity, which is √ −τ w = Π H/ρ f , with τ 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 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 = 11 250), Pr = 1 and St = 1. For each Ri * , we set T t such that the initial heat flux in the ice, −T t /2, is almost equal to the heat flux in the fluid when there is no melting. As a result, the fluid-solid interface position does not move significantly in time (at least initially) and we can maximize numerical resolution around the interface with a fixed grid. For reference, the Rayleigh number for the unstable stratification case (Ri * = −40) is Ra = 4.5 × 10 5 , which is above the instability onset for Rayleigh-Bénard convection rolls in the streamwise direction (Ra ≈ 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 α, which indeed can be either positive or negative depending on the fluid state. We discuss the geophysical relevance of our choice of parameters, including α, in § 4.
We solve (2.5) using the open-source pseudo-spectral direct numerical simulation (DNS) code DEDALUS . 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 Δz + = 1/Re * ≈ 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 φ = 1 in the fluid to φ = 0 in the solid occurs and Nu III b 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; ξ III , ξ − and ξ + 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 I D and C III D 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 × 10 5 . over a very thin diffuse interface of thickness ≈0.007. Simulation parameters and output variables are provided in table 1.

Variables of interest
We define the friction velocity, the bulk velocity and the Nusselt number as respectively (cf. details in appendix B), where the overbar denotes the horizontal average and · ≡ V dV/V denotes the volume average, such that φ is the mean fluid fraction and u b is the bulk velocity of the fluid phase. In (2.9a-c), τ d , τ ν and τ w are the linear damping, viscous and Reynolds shear stresses, and q = wT − ∂ z T is the heat flux. At statistical steady state, the full shear stressτ = (τ d + τ ν + τ 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τ is a linear function of z, u * can be estimated from the full shear stress as √ −τ 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 φ ≤ 1, such that u * is truly the friction velocity at the mean interface position (cf. (2.9a-c)). We denote by ξ the fluid-solid interface position, where (2.10) such thatξ = φ (note that one could alternatively define ξ such that it satisfies φ(z = ξ) = 0.5 or T(z = ξ) = 0), and we denote the melt rate byṁ = ∂ t ξ . 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 2 b /2, i.e. 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 χ(t, x, y) = ξ(t, x, y) − z.

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.

Simulation stages and global flow variables
We show in figure 2 the friction velocity u * , the bulk velocity u b and the Nusselt number Nu for stable (top figure), neutral (middle figure) and unstable (bottom figure) stratification. Each simulation is broken down into three main stages, which are highlighted by different colours in figure 2 (note that we do not discuss the results shown in grey, which correspond to the spin-up of the fluid phase without buoyancy effects). The first stage of interest (stage I for t Ib * ≤ t ≤ t Ic * ) is shown in blue and corresponds to the spin-up of the fluid phase with buoyancy effects turned on. Importantly, stage I neglects the solid phase, which is substituted with a simple isothermal no-slip boundary, for computational expediency. The second key stage (stage II for t Ic * < t ≤ 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 ≤ z ≤ 1.5 and the phase field is prescribed as φ = 0.5 1 − tanh 2(z − 1)/δ , where δ 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 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 with Δ * = 50, and are reported with q s in table 1. Note that Re b 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 * ≤ t Ib * , such that the results of figure 2 are exactly the same for all three simulations until t * = t Ib * . Upon turning on buoyancy, i.e. for t * ≥ t Ib * (blue), the Nusselt number and bulk velocity deviate from the neutral case (middle figure), but with opposite behaviours: Nu decreases while u b increases with stabilizing buoyancy effects (top figure), and Nu increases while u b decreases with destabilizing buoyancy effects (bottom figure). The friction velocity, on the other hand, remains close to u * /Pr ≈ 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 I = 3.02, 4.58 and 6.74 for Ri * = 40, 0 and −40, respectively (cf. table 1). It is worth noting that while Nu remains the same between stage I and stage II (in a time-average sense), u * and u b show some variations as a result of turning on volume penalization and adding a solid phase. The large dip of u * at t * ≈ t Ic * is merely the result of a sudden deceleration of the mean flow close to the interface, due to the addition of linear damping, which is transient, as can be seen from the rapid return of u * to its statistically steady-state value of u * ≈ 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 Ib * , i.e. when the stratification becomes stable (respectively unstable). The decrease (respectively increase) of C D results from an increase (respectively decrease) of the potential energy barrier in stirring the mean shear and bringing momentum upward with increasing stable (respectively unstable) stratification and is in agreement with previous studies (García-Villalba & del Álamo 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).

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 ( ) and variables averaged in the x direction minus the horizontal mean by a prime ( ), such that e.g. ξ =ξ −ξ represents the spanwise variations of the streamwise-averaged topography around the horizontal mean.
We show snapshots of the two-dimensional fluid-solid interface ξ 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 ξ 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, ξ + (maximum of ξ ), and the steady-state amplitude of the deepest keels, ξ − (minus the minimum of ξ ), 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 δ ν = 1/Re * for stable, neutral and unstable stratification, respectively (note that δ ν 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δ ν , 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  with unstable stratification. For stable stratification (figure 4d), the separation of scales between the topography lifetime (about 10 friction time units) and the diffusion time scale across the solid layer (about 100 friction time units) suggests that the interface evolution is purely driven by the flow dynamics. For neutral stratification, figure 4(e) shows that channels and keels can drift, merge, split, decay and spontaneously appear over time scales of tens to hundreds of friction time units, highlighting a possible interplay between interface evolution and the fixed-temperature condition at the top solid boundary. For unstable stratification (figure 4f ), the channels and keels become time invariant and their amplitudes saturate because of the top solid boundary condition, which plays a key role in the interface evolution as discussed in the next sections.

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   in the case of unstable stratification (figure 5c), −∂ zT 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 Ri * . In the case of stable stratification (figure 5a), buoyancy effects inhibit the generation of large topographical features, such that channels and keels have small amplitudes and do not feed back onto the flow (which is further discussed in § 3.4). In the case of neutral stratification, buoyancy is turned off, such that the solid boundary deforms more and can affect the flow dynamics. Figure 5(b) shows that the heat flux in the fluid close to the boundary is usually larger where there are channels (e.g. y ≈ 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 * ≈ 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 * ≥ 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 ofq (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π) = 2, which is close to the critical wavenumber L y k c /(2π) = 3.11 of convection instability (λ 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 ξF igure 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ξ 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 − χ, where χ 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. 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ξ 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ξ 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 χ ) 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 χ < 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 = ξ − 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 911 A44-20  coming from the top boundary. Also, for the vertical velocity at z = ξ − 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 , such that Ri * starts from ≈−40 at t * = t III * and reaches ≈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 (ṽ,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 * ≈ −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 * ≈ 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 * ≈ 0.7 and Ri * ≈ 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.

Melting
In this section we discuss the evolution of the mean interface positionξ with time and the statistics of meltingṁ = ∂ t ξ at the statistical steady state. We first show the evolution ofξ 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ξ with time for an unstable than for a stable stratification, since the heat flux in the fluid is larger when buoyancy forces are destabilizing rather than restoring. However, as indicated in § 2, we have imposed a conductive heat flux in the solid (q s ) slightly less than the mean heat flux through the fluid (Nu I ) at the beginning of stage III, i.e. when we turn on melting, such that the leading-order melt rate is not controlled by the Nusselt number of the fluid-only simulations but by the difference Nu I − q s . This difference is 0.19, 0.1 and 0.05 for stable, neutral and unstable stratification (cf. table 1), such that the initial increase ofξ is faster for stable than for unstable stratification.  Figure 11 shows thatξ saturates over time. This happens because the conductive heat flux in the ice increases asξ 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ξ becomes where h 0 = 1/2 is the initial ice thickness, q f is the heat flux in the fluid and we recall thatξ = ξ 0 = 1 at t = t II * . For simplicity, we take q f to be a constant diagnosed from the simulations. The results of (3.4) for q f = Nu I , i.e. obtained when setting q f to the average heat flux before melting is turned on, are shown by the dotted lines in figure 11(a). The overlap between the reduced model and DNS results at early times is good for unstable stratification (as expected) but is poor for stable and neutral stratification. The disagreement with q f = Nu I at early times arises because the temperature in the solid is slightly above 0 because of volume penalization, such that there is some artificially large melting at the beginning of stage III (cf. appendix A). At later times, the DNS results and the model results shown by the dotted lines diverge because the heat flux Nu increases rapidly once melting is turned on, as can be seen in figure 11(b). For unstable stratification, the agreement with q f = Nu I is relatively good until t * ≈ t II * + 50 (cf. red dotted line), i.e. right until the Rayleigh-Bénard rolls are energized and a large-scale topography emerges (cf. § 3.3).
In order to account for the increase in heat flux through the fluid enabled by melting and the generation of topography, we show with dashed lines in figure 11(a) the result of (3.4) with q f = Nu III , which is the heat flux at the statistical steady state with melting turned on. For stable and neutral stratification, there is good agreement between the model results and the mean interface position at late times. For unstable stratification, however, (3.4) with q f = Nu III overestimates the final value ofξ − 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ξ 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ṁ are shown in figure 12(b). The temporal and spatial average,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 (ṁ < 0), but occasionally melts rapidly (ṁ > 0). We remark that the asymmetry of the melt rate p.d.f.s is not due to the asymmetry of the interface p.d.f.s since the p.d.f.s of melt rates inside channels (upper triangle) and along keels (lower triangles) are similar, but is instead a generic feature of melting by a turbulent fluid. Indeed, while the turbulent flow can drive rapid melting independently of what happens in the solid, freezing necessarily involves slow diffusive processes in the solid. Additionally, the near-wall dynamics, which features coherent structures such as streamwise streaks and vortices, is itself asymmetric, as can be seen from the p.d.f.s of the temperature gradient in figure 9. Thus, it is not surprising that the melt rate p.d.f.s are asymmetric. While beyond our goal, it would be worthwhile in the future to try to identify flow features controlling the shape of the melt rate p.d.f.s.

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. L = 3 × 10 5 J kg −1 and c p = 4 × 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 ν = κ = 10 −6 m 2 s −1 , i.e. considering an anomalously high thermal diffusivity such that Pr = 1, then our simulation runs correspond to about 10 hours in real time, the bulk velocity is 2 cm s −1 and the thermal expansion coefficient is 5 × 10 −7 K −1 in absolute value, which is quite small for water. These calculations highlight that there is clearly a significant gap between the control parameters in our simulations and in nature. Our set of experiments is also limited to three different runs, such that we cannot offer a quantitative prediction of what would be observed in the environment. Nevertheless, Rayleigh-Bénard rolls have been observed for a wide range of Rayleigh numbers and Reynolds numbers in mixed-convection simulations with unstable stratification (Pirozzoli et al. 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 × 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 2 b ≈ 0.1 (defined as positive for an unstable stratification) and large Monin-Obukhov length L MO ≈ 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 ≥ 0.1, and to grow in size as Ri b increases. In the limit Ri b → ∞, channels and keels may be expected to eventually disappear. Indeed, as Ri b → ∞, 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 * ≥ 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 * ≥ 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 λ usually satisfies Re λ * = Re * λ/H ≥ O(1000-10 000) and have always reported a scallop wavelength smaller than the fluid depth, i.e. λ < H (Blumberg & Curl 1974;Thomas 1979;Claudin et al. 2017). Considering the upper limit λ = H means that scallops are predicted to emerge for Re * ≥ O(1000-10 000) in water, which is at least one order of magnitude higher than what we selected for our study and difficult to achieve numerically. Note, however, that the minimum Re * leading to scallop formation may be different for a fluid with control parameters Pr = 1 and St = 1 (as is the case in this work) instead of Pr ≈ 10 and St ≈ 75, as is the case for water, and may also vary with the stratification strength.

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 3 with a spectral code) would still be required. It is noteworthy that simulations of a pure fluid at a fixed wavy boundary would already be useful in helping to verify or refine the most recent theoretical predictions on scallop formation and saturation (Claudin et al. 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 * 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 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. A = 6/(5St), B = (16/δ 2 ) × 6/(5St) and C = 1, with δ chosen such that it is equal to two times the local grid size at z = 1 (initial interface position). Also, Γ = (δ/2.648228) 2 and we require time steps to be always smaller than Γ /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 δ is the thickness of the diffuse phase-field interface over which φ transitions from 1 in the fluid to 0 in the solid (see figure 1c). Here δ 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 δ ν = 1/Re * . Here, we chose to have δ ≈ δ ν in all simulations, i.e. δ is equal to 1 wall unit Δz + = δ, 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 δ ν /2. The damping time scale Γ in (2.5) is set to Γ = (δ/2.648228) 2 in order to cancel first-order errors in the phase-field model . 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 + ≥ 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 + ≈ 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 + | < δ (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 + < δ 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 -10 0 0 10 0 10 1 10 2 -10 0 0 10 0 10 1 10 2 -10 0 0 10 0 10 1 10 2 Figure 14. Wall-normalized velocity U + = Pr u/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 * ∈ [t Ib * , t Ic * ]) whereas orange symbols denote results obtained with volume penalization (i.e. for t * ∈ [t Ic * , t II * ]). The plus symbols show the temperature profiles obtained with volume penalization but shifted to the right, i.e. into the fluid, by replacing z + with z + + 2.5, z + + 1.16 and z + + 0.78 for stable, neutral and unstable stratification. The overbars denote horizontal averaging and time averaging over 20 friction time units. The vertical dotted lines show z + = −δ, 0, δ, where δ is the diffuse interface thickness.
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 + ≈ 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).
whereτ andq are the anomalous stress and heat flux due to damping of the advective terms in (2.5), i.e.
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 τ w , the viscous stress τ ν , the linear damping stress τ d and the Reynolds stress plus the anomalous stress τ w +τ for stable, neutral and unstable stratification, respectively. Importantly, τ w and τ w +τ 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 τ w + τ ν + τ d +τ , i.e. the total stress (thin dashed lines) that includes the anomalous stress. We show in figure 15(g-i)q (solid lines) and q +q (thin dashed lines). For stable and neutral stratification (top and middle rows),q and q +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,q is not perfectly constant, but deviates from q +q and peaks at z ≈ 1.15, which is roughly the height of the channels. The relative discrepancy betweenq and q +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. (u · ∇)(φ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 ≤ t Ic * (including stage I), we solve (2.5a), (2.5b) and (2.5d) with φ ≡ 1 and a no-slip isothermal top-boundary condition, i.e. u = 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 φ prescribed, i.e. not varying in time (cf. the main text). In stage III we solve (2.5) with all variables freely evolving and we use the same spectral resolution as in stage II.
Our simulations until t = t Ic * can be broken down into three sub-stages. In stage Ia, i.e. for t * < t Ia * (cf. light grey colour in figure 2), we run a low-resolution (128 Fourier modes in the x and y directions and 32 Chebyshev modes in the z direction) spin-up simulation of an initially laminar flow superposed with three-dimensional velocity perturbations and no buoyancy effects (Ri * = 0). In stage Ib, i.e. for t Ia * ≤ t * < t Ib * , we increase the resolution 911 A44-32 -20 000 20 000 0 -20 000 20 000 0 Figure 15. (a-c) Reynolds stress τ w (dotted lines), viscous stress τ ν (dashed lines), linear damping stress τ d (solid lines) and Reynolds stress plus anomalous stress τ w +τ (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. τ w + τ ν + τ d (solid lines), and total stress plus the anomalous stress, i.e. τ w + τ ν + τ d +τ (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 +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. and the second-order solution is Thus, the second-order-accurate formula for the mean heat flux at the top of the ice reads as 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 = 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.