2.1. The one-phase Stefan problem

We consider a 1-D Stefan problem in which an ice layer of thickness $\xi$ is subjected to an arbitrary constant thermal forcing, $q_w$ , at the bottom boundary (corresponding to the moving ice–water interface) and an imposed temperature on the top boundary. The ice layer extends from $z=0$ to $z = \xi (t)$ , which represents the instantaneous position of the ice–water interface. Any point where $z\gt \xi (t)$ falls inside the water layer. The temperature, $\theta = \theta (z,t)$ , is set equal to a constant value below the freezing point, $\theta _c \lt 0$ , at $z=0$ and to the melting temperature, $\theta _m=0$ , at $z = \xi (t)$ . Moreover, $ q_w$ is defined such that it is positive when transferred from ice to water: hence, $ q_w$ is always negative during melting. The dimensionless governing equations (see Appendix A for details on their non-dimensionalization) read as

(2.1) \begin{equation} \left \{\!\! \begin{array}{l} \displaystyle {{\textit{Pe}}\,{St}\dfrac {{\textrm {d}} \xi }{{\textrm {d}} t} = \left .\dfrac {\partial \theta }{\partial z}\right |_{z=\xi } + q_w} ,\\[2ex] \displaystyle {\dfrac {\partial \theta }{\partial t}=\dfrac {1}{{\textit{Pe}}}\dfrac {\partial ^2 \theta }{\partial z^2}} , \end{array} \right . \; \end{equation}

with the following boundary conditions

(2.2) \begin{align} \theta (z=0,t)&=\theta _c , \end{align}

(2.3) \begin{align} \theta (z=\xi (t),t)&=0 , \end{align}

(2.4) \begin{align} \theta (z,t=0)&=\theta _c(1-z/\xi_0) , \end{align}

(2.5) \begin{align} \xi (0)&=\xi _0 \, . \end{align}

The dimensionless groups in (2.1) are defined as follows

(2.6) \begin{equation} {St}=\frac {\tilde {\mathcal{L}}}{ \tilde {c}_p \Delta \tilde {\theta }} \; , \qquad {\textit{Pe}}=\frac {\tilde {\rho } \tilde {c}_p \tilde {u}_{\textit{ref}}\, \tilde {h}_{\textit{ref}}}{\tilde {\lambda }} \, ; \end{equation}

where the notation $\tilde {\bullet }$ indicates that the quantity $\bullet$ is dimensional; $\tilde {\mathcal{L}}$ is the specific latent heat, $\tilde {\rho }$ the density, $\tilde {c}_p$ the specific heat at constant pressure, $\tilde {\lambda }$ the thermal conductivity, $\Delta \tilde {\theta }$ the reference temperature difference across the ice layer (defined here as $\Delta \tilde {\theta } = \tilde {\theta }_m-\tilde {\theta }_c$ ), $\tilde {u}_{\textit{ref}}$ the reference velocity and $\tilde {h}_{\textit{ref}}$ the reference length of the problem. In a real case scenario, heat diffusion is significantly faster than melting. This means that the temperature distribution within the ice layer can quickly adapt to the motion of the interface and can always be considered in equilibrium state. It is therefore reasonable to assume that the time scale of heat diffusion, $\tau _{\mathcal{D}}$ , is much smaller that the time scale of melting, $\tau _{\mathcal{M}}$ . These two time scales can be defined as

(2.7) \begin{equation} \left \{\!\! \begin{array}{l} \displaystyle {\left |\dfrac {{\textrm {d}} \xi }{{\textrm {d}} t}\right | =\dfrac {\xi }{\tau _{\mathcal{M}}}} , \\[2ex] \displaystyle {\left |\dfrac {\partial \theta }{\partial t}\right |=\dfrac {\Delta \theta }{\tau _{\mathcal{D}}}} , \\ \end{array} \right . \; \end{equation}

where $\varDelta$ is the characteristic temperature difference inside the ice layer (note that $\Delta \theta = 1$ by definition). The condition $\tau _{\mathcal{M}} \gg \tau _{\mathcal{D}}$ , in combination with the system (2.1), yields

(2.8) \begin{equation} {St} \gg \dfrac {\Delta \theta \left |\dfrac {\partial \theta }{\partial z}+ q_w\right |}{\xi \left |\dfrac {\partial ^2 \theta }{\partial z^2}\right |}\approx \Delta \theta \dfrac {\Delta \theta / \xi + q_w}{\Delta \theta / \xi } \; \; \implies \; \; {St} \gg \left | 1 + q_w \xi \right | . \end{equation}

When the condition (2.8) is satisfied, the following analytical solution for the problem (2.1) exists

(2.9) \begin{equation} \left \{\!\! \begin{array}{l} \displaystyle {\xi (t) = \dfrac {\theta _c}{ q_w}}\dfrac {S-1}{S} = \xi _{\textit{eq}} \dfrac {S-1}{S} , \\[2ex] \displaystyle {\dfrac {1}{S} + \ln {\left |S\right |} = \dfrac {1}{S_0} + \ln {\left |S_0\right |} - \dfrac { q_w^2}{\theta _c}\dfrac {t}{{\textit{Pe}}\,{St}}} , \\[2ex] \displaystyle {S_0 = \dfrac {1}{1- \xi _0/\xi _{\textit{eq}}}} , \\[2ex] \theta (z,t) = \theta _c\left (1-z/\xi (t)\right ) , \end{array} \right . \; \end{equation}

where $\xi _{\textit{eq}}=\theta _c/ q_w$ is the ice thickness at equilibrium. The reader is referred to Appendix B for a detailed derivation of this solution. The second equation in (2.9) can be used to obtain an estimate of the time required to reach the steady state, $\mathcal{T}_{\textit {steady}}$ . In principle, this equation predicts an infinite time since $S \rightarrow \infty$ (and $\ln {|S|} \rightarrow \infty$ ) for $\xi = \xi _{\textit{eq}}$ . Therefore, we define $\mathcal{T}_{\textit {steady}}$ as the time required to melt $99 \, \%$ of all the ice that can melt once the equilibrium condition $\xi = \xi _{\textit{eq}}$ is reached. When this condition is met, $\xi = \xi _0 - 0.99(\xi _0-\xi _{\textit{eq}})$ , which yields $S^{-1} = 0.01(\xi _{\textit{eq}}-\xi _0)/\xi _{\textit{eq}}$ and, in turn,

(2.10) \begin{equation} \mathcal{T}_{\textit {steady}} = \dfrac {\theta _c\,{\textit{Pe}}\,{St}}{ q_w^2}\left [ 0.99(1- q_w\xi _0/\theta _c) + \ln {(0.01)} \right ] \! . \end{equation}

2.2. The two-phase Stefan problem

To complete the model and account for the heat transfer balance in the water layer, the heat flux $q_{S,w}$ that brings the meltwater to local thermal equilibrium with the bulk water must be included. This contribution is especially important at low St, when it is comparable to the latent heat $q_{L}$ and cannot be neglected. In fact, if the assumption $q_{S,w}\approx 0$ is made (which is analogous to assuming a constant heat flux $q_w$ applied directly to the interface, regardless of the velocity with which it evolves), the model predicts an instantaneous evolution of $\xi$ to the equilibrium position in the limit ${St} \to 0$ . In other words, the model predicts an infinite velocity of the ice–water interface at vanishing Stefan numbers. In reality, for small but finite values of $ {St}$ , the interface velocity is limited by thermal diffusivity, which is responsible for the presence of $q_{S,w}$ . In the case of ice melted by a turbulent flow of water, $q_{S,w}$ can be quantified as follows (in physical units)

(2.11) \begin{equation} \tilde {q}_{S,w}=\tilde {q}_b-\tilde {q}_w = \big(\tilde {\theta }_{b}-\tilde {\theta }_m\big)\tilde {\rho } \tilde {c}_p\dfrac {{\textrm {d}} \tilde {\xi }}{{\textrm {d}} \tilde {t}} , \end{equation}

with $\tilde {\theta }_{b}$ the bulk temperature of water. Equation (2.11) shows that $\tilde {q}_{S,w}$ can be interpreted as the heat flux required to bring the meltwater layer from the melting temperature $\tilde {\theta }_m$ to the bulk temperature $\tilde {\theta }_{b}$ . In dimensionless form, (2.11) becomes

(2.12) \begin{equation} q_{S,w} = q_b- q_w = \underbrace { \left [ \dfrac {\big(\tilde {\theta }_{bulk}-\tilde {\theta }_m\big)\tilde {\rho } \tilde {c}_p}{\tilde {\rho } \tilde {c}_p \Delta \tilde {\theta } } \right ] }_{\mathcal{Q}} \underbrace { \left ( \dfrac {\tilde {\rho } \tilde {c}_p \tilde {u}_\tau }{\tilde {\lambda }} \right ) }_{{\textit{Pe}}} \dfrac {{\textrm {d}} \xi }{{\textrm {d}} t} = \mathcal{Q} \; {\textit{Pe}} \; \dfrac {{\textrm {d}} \xi }{{\textrm {d}} t} , \end{equation}

where $\tilde {\rho } \tilde {c}_p \Delta \tilde {\theta }$ is the reference sensible heat. Substituting (2.12) into (2.1) and assuming constant $q_b$ yields a problem equal to that solved in the previous section, provided that St is now replaced by $\mathcal{Q} + {St}$ and $ q_w$ by $ q_b$ . Therefore, the analytical solution for the evolution of the ice layer thickness, $\xi$ , is

(2.13) \begin{equation} \left \{\!\! \begin{array}{l} \displaystyle {\xi (t) = \xi _{\textit{eq}}}\dfrac {S-1}{S} ,\\[2ex] \displaystyle {\dfrac {1}{S} + \ln {\left |S\right |} = 1-\dfrac {\xi _0}{\xi _{\textit{eq}}}-\ln {\left |1-\dfrac {\xi _0}{\xi _{\textit{eq}}}\right |} - \dfrac { q_b^2}{\theta _c}t^*} , \end{array} \right . \end{equation}

where $t^*=t/ [{\textit{Pe}} ({St}+\mathcal{Q} ) ]$ is the rescaled time. The expression for the time required to reach the steady state, stemming from (2.10), becomes

(2.14) \begin{equation} \mathcal{T}_{\textit {steady}} = \dfrac {\theta _c{\textit{Pe}}\left ({St}+\mathcal{Q}\right )}{ q_b^2}\left [ 0.99(1- q_b\xi _0/\theta _c) + \ln {0.01} \right ]\!. \end{equation}

The accuracy of the model (2.13) has been verified by comparing its predictions with the numerical solution of the full set of (2.1). Details of the model assessment are provided in Appendix B.1.

3.1. Governing equations

The system of dimensionless governing equations reads as

(3.1) \begin{equation} \left \{\!\! \begin{array}{l} \displaystyle {\frac {\partial \phi }{\partial t}=\frac {6}{5}\frac {1}{\mathcal{C}\,{\textit{Pe}}\,{St}}\left [\nabla^{2}\phi -\frac {1-\phi }{\epsilon ^2}\phi \left (1-2\phi +{\mathcal{C}} { \theta }\right )\right ]} , \\[2ex] \displaystyle {\frac {\partial { \theta }}{\partial t}=\frac {1}{{\textit{Pe}}}\nabla^{2} { \theta }-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\theta +{St}\frac {\partial \phi }{\partial t}} , \\[2ex] \displaystyle {\frac {\partial \boldsymbol{u}}{\partial t} =\frac {1}{{\textit{Re}}_{\tau }}{\nabla} ^2\boldsymbol{u}-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}-\boldsymbol{\nabla }p-\frac {\phi ^2}{\eta _s}\boldsymbol{u}} , \\[2ex] \displaystyle {\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0} \, . \end{array} \right . \; \end{equation}

The first equation in this system is the modified Allen–Cahn (AC) equation, used to track the evolution of the local ice volume fraction, $\phi$ . The second equation is the EN, used to track the evolution of the temperature field, $\theta$ . The last two equations are the Navier–Stokes and continuity equations, used to track the evolution of the velocity field, $\boldsymbol{u}$ . The dimensionless groups that appear in (3.1) are the shear Reynolds number, the Péclet number and the Stefan number, defined as

(3.2) \begin{equation} {\textit{Re}}_{\tau }=\frac {\tilde {u}_\tau \tilde {h}}{\tilde {\nu }} \; , \qquad {\textit{Pe}}={\textit{Re}}_{\tau }\,{{\textit{Pr}}} =\frac {\tilde {\rho } \tilde {c}_p \tilde {u}_{\tau } \tilde {h}}{\tilde {\lambda }} \; , \qquad {St}=\frac {\tilde {\mathcal{L}}}{\tilde {c}_p \Delta \tilde {\theta }} \; , \end{equation}

with $\tilde {\nu }$ the kinematic viscosity, $\tilde {\rho }$ the density, $\tilde {c}_p$ the specific heat and $\tilde {\lambda }$ the thermal conductivity. For simplicity, we consider constant thermophysical properties across the two phases. The definitions of St and ${\textit{Pe}}$ given in (3.2) are analogous to those provided in § 2, with $\tilde {u}_{\tau }$ replacing $\tilde {u}_{\textit{ref}}$ and $\tilde {h}$ replacing $\tilde {h}_{\textit{ref}}$ . This yields ${\textit{Pe}}=Re_\tau Pr$ with ${{\textit{Pr}}}=(\tilde {\rho }\tilde {c}_p\tilde {\nu })/\tilde {\lambda }$ the Prandtl number. The flow of water is driven by a mean pressure gradient imposed in the streamwise direction. Since the average height of the water layer, $\tilde {h}_w(t)$ , increases over time due to melting, the effective Reynolds number of the flow, ${\textit{Re}}_{\textit{eff}} = \tilde {u}_\tau \tilde {h}_w/\tilde {\nu }$ , changes accordingly and ${\textit{Re}}_{\tau }$ corresponds to the maximum value that ${\textit{Re}}_{\textit{eff}}$ can reach (when $\tilde {h}_w = \tilde {h}$ ). In the following, we indicate the initial value of the Reynolds number as ${\textit{Re}}_{\tau ,0}$ and its final value as ${\textit{Re}}_{\tau ,f}$ , such that ${{\textit{Re}}_{\tau ,0}} \leqslant {\textit{Re}}_{\textit{eff}} \leqslant {\textit{Re}}_{\tau ,f}$ . In the type of problem we decided to investigate, temperature is treated as a passive scalar and therefore buoyancy effects are negligible (Zonta & Soldati Reference Zonta and Soldati2018).

The modified AC equation is used to implicitly describe the morphodynamic changes of the ice–water interface and is derived from a formulation of the phase-field method specifically designed for problems involving melting and solidification (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020). The equation includes a source term on the right-hand side that depends on the local temperature $\theta$ , such that if the local temperature is below the freezing point, then ice is generated, while water is generated and the ice volume fraction is reduced if $\theta$ is above the melting point. The phase field method ensures a smooth transition between the ice phase (where $\phi =1$ ) and the water phase (where $\phi =0$ ) across a thin separation layer centred around the ice–water interface (where $\phi =0.5$ ). The thickness of this layer is dictated by the phase field parameter, $\epsilon$ , which is imposed by computational requirements. In our simulations, $\epsilon =0.005$ , a value that corresponds to a layer approximately $0.04 \tilde {h}$ thick. The parameter $\mathcal{C}$ appearing in the AC equation is the mobility coefficient of the phase-field method. From a physical point of view, this parameter is associated with the energy of the transition layer that separates the two phases. This energy (and, therefore, the value of $\mathcal{C}$ ) affects the dependency of the temperature at the interface on the local curvature in conditions of thermodynamic equilibrium (Gibbs–Thompson effect). Following Yang et al. (Reference Yang, Howland, Liu, Verzicco and Lohse2023b ) and Perissutti et al. (Reference Perissutti, Marchioli and Soldati2024), we neglect the Gibbs–Thompson effect, since the curvatures of the ice–water interface are never large enough to make the interfacial temperature significantly different from the melting point. Numerically, this is achieved by setting ${\mathcal{C}}=10$ (Perissutti et al. Reference Perissutti, Marchioli and Soldati2024). This value does not affect significantly the numerical stiffness of the problem, which increases with $\mathcal{C}$ . In such conditions, the phase field formulation enforces the temperature to always be at the melting/freezing point ( $\theta =0$ , as it follows from the definition of dimensionless temperature) at the ice–water interface. In the EN, a source term proportional to the time derivative of $\phi$ and to St is introduced to account for the contribution of the latent heat during melting and freezing. Following previous studies (Favier, Purseed & Duchemin Reference Favier, Purseed and Duchemin2019; Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020; Couston et al. Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2021; Yang et al. Reference Yang, Chong, Liu, Verzicco and Lohse2022; Perissutti et al. Reference Perissutti, Marchioli and Soldati2024), a volume-penalization immersed boundary method (IBM) is used in combination with the Navier–Stokes equations. In this method, a source term is added in the momentum equations to force an exponential decay of the fluid velocity where the ice volume fraction is not zero. This decay is ruled by the time constant $\eta _s$ , chosen as small as possible ( $7\times 10^{-4}$ in this study) to minimize the error between the exact solution of the Navier–Stokes equations and the IBM solution, which scales as $\sqrt {\eta _s}$ when compared with the sharp interface formulation (Favier et al. Reference Favier, Purseed and Duchemin2019), and prevent any velocity in the ice layer.

3.2. Numerical discretization and boundary conditions

The Navier–Stokes and continuity equations in the system (3.1) are rewritten in the wall-normal velocity–vorticity formulation (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2007) and solved, together with the AC and EN equations, using a pseudospectral method that relies on discrete Fourier transforms along the homogeneous directions ( $x$ and $y$ ) and discrete Chebyshev transforms along the wall-normal direction ( $z$ ). For this reason, the computational grid (consisting of $N_x\times N_y \times N_z$ points) is uniform in $x$ and $y$ , while Chebyshev–Gauss–Lobatto points are used in the $z$ direction. The equations are advanced in time using an implicit–explicit scheme: the nonlinear terms are discretized using an Adams–Bashforth scheme, while for the linear terms, an implicit Euler scheme (for AC and EN) or a Crank–Nicolson scheme (for the wall-normal velocity–vorticity equations) is employed. More details about the numerical scheme can be found in Roccon (Reference Roccon2024) and Roccon, Soligo & Soldati (Reference Roccon, Soligo and Soldati2025). The solver is parallelized using a message passing interface (MPI) and leverages on a two-dimensional domain decomposition strategy. The code is also accelerated exploiting graphics processing units (GPUs) through OpenACC directives and CUDA Fortran instructions. The Nvidia cuFFT libraries are used to execute the Fourier/Chebyshev transforms.

As far as the boundary conditions are concerned, for the Navier–Stokes equation, a no-slip condition is enforced at the top boundary capping the ice layer, while a free shear condition is applied at the bottom boundary of the water layer,

(3.3) \begin{equation} { \left .\boldsymbol{i}\boldsymbol{\cdot }\frac {\partial \boldsymbol{u}}{\partial z}\right |_{z = 0}=0 , \qquad \left .\boldsymbol{j}\boldsymbol{\cdot }\frac {\partial \boldsymbol{u}}{\partial z}\right |_{z = 0}=0 , \qquad \boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{u}(z=0)=0 , \qquad \boldsymbol{u} (z=1)=\boldsymbol{0}\, . } \end{equation}

Note that the velocity is zero throughout the ice layer since a no-slip condition is enforced at the ice–water interface by the IBM. In the EN, constant temperatures are imposed at both boundaries: above the melting point at the bottom boundary, $\theta (z=1)=1$ , and below the freezing point at the top one, $\theta (z=0)=-1$ . The choice imposing a symmetric temperature difference between the top and bottom boundaries was made to be fully consistent with our previous work (Perissutti et al. Reference Perissutti, Marchioli and Soldati2024) and ensures that the ice–water interface remains inside the domain. We remark here that imposing an asymmetric temperature difference would only change the vertical equilibrium position of the ice interface inside the domain: no significant changes would be produced in the ice morphology. For the AC equation, no-flux conditions are applied at both boundaries,

(3.4) \begin{equation} \left .\frac {\partial \phi }{\partial z}\right |_{z= 0}=0 , \qquad \left .\frac {\partial \phi }{\partial z}\right |_{z= 1}=0 \, . \end{equation}

These conditions prevent any unphysical flux of the ice phase through the domain boundaries and ensure that the changes in the volume fraction $\phi$ are only due to melting and freezing. Along the homogeneous directions ( $x$ and $y$ ), periodic boundary conditions are applied to all flow variables.

3.3. Simulation set-up

All simulations, summarized in table 1, start from the same initial condition, which consists of an initially flat ice layer in the top part of the domain, from $z=0.75$ to $z = 1$ . The initial volume fraction is

(3.5) \begin{equation} \phi _0(x,y,z) = \dfrac {1}{2}\tanh \left (\dfrac {z-0.75}{2\epsilon }\right )+\dfrac {1}{2} , \end{equation}

corresponding to the phase field equilibrium profile in the absence of melting and freezing. The interface position associated with (3.5) ensures that the interface is not too close to the channel centre (where the grid is coarser because of the Chebyshev nodes distribution), but the ice layer is still thick enough to properly catch the transient of the melting process. Below the flat ice layer, a fully developed turbulent flow is initialized by running a preliminary simulation in which ice melting is artificially suppressed, and the ice–water interface is kept flat and its thickness unchanged. This preliminary simulation is carried out until the statistically steady state is reached for the flow field. The steady-state flow is then used as the initial condition for the simulations with melting.

Table 1.

Overview of the simulation parameters for each simulation (S1 to S6): Reynolds (initial and final), Stefan and Prandtl numbers, number of grid points ( $N_x$ , $N_y$ , $N_z$ ), the ice morphology regime (subcritical or supercritical), the predicted time required for the average ice thickness to reach its steady-state value, $\mathcal{T}_{\textit {steady}}$ , and the estimated computational cost, normalized by the cost of the cheapest simulation, S1. One computational cost unit is roughly 50 GPU hours on a petascale Tier-0 supercomputer.

To conduct the study, two Reynolds numbers were simulated: one in the absence of the instability of the ice–water interface that leads to the formation of ice ripples (referred to as subcritical regime hereinafter, corresponding to a shear Reynolds number at the beginning of the simulation ${\textit{Re}}_{\tau ,0}=170$ ) and one in which ice ripples form as a result of the interaction between ice streamwise advection and cross-stream diffusion of heat and momentum (referred to as supercritical regime hereinafter, corresponding to ${\textit{Re}}_{\tau ,0}=636$ ). These values were chosen considering that the characteristic wavelength $\lambda ^+_{x,cr}$ of the unstable disturbances in the flow should be larger than 2100 in wall units, with the fastest growth rate being achieved at $\lambda ^+_{x,cr} \simeq 3500$ (Ashton & Kennedy Reference Ashton and Kennedy1972; Hsu et al. Reference Hsu, Locher and Kennedy1979). The corresponding (estimated) critical Reynolds numbers above which ice ripples are expected to form, are: ${\textit{Re}}_\tau ^{cr} = 251$ to have unstable modes, and ${\textit{Re}}_\tau ^{cr} = 418$ to have the fastest growing mode (Perissutti et al. Reference Perissutti, Marchioli and Soldati2024). It is worth mentioning that, at present, it is not possible to quantitatively predict a critical Reynolds number to discriminate between the two regimes based solely on the physical parameters of the system, primarily due to the difficulty of modelling with accuracy the phase shift between the heat flux and the ice–water interface (Perissutti et al. Reference Perissutti, Marchioli and Soldati2024). The phase shift is typically taken into account using empirical corrections (Thorsness, Morrisroe & Hanratty Reference Thorsness, Morrisroe and Hanratty1978), although recent studies have attempted to address this limitation (Chedevergne et al. Reference Chedevergne, Stuck, Olazabal-Loumé and Couzi2023).

Table 1 shows ${\textit{Re}}_{\tau ,0}$ together with the effective Reynolds number at the end of the simulation ( ${\textit{Re}}_{\tau ,f}$ ). For each value of the Reynolds number, we simulated three different Stefan numbers to evaluate the effect of the latent heat contribution. We already anticipated in § 2 that the main effect of St is to change the time required to reach the steady state, $\mathcal{T}_{\textit {steady}}$ and that, at large St, $\mathcal{T}_{\textit {steady}}$ (and therefore the computational cost) becomes significantly higher. For this reason, table 1 also shows the value of $\mathcal{T}_{\textit {steady}}$ for each simulation as obtained from (2.14). Note that $\mathcal{T}_{\textit {steady}}$ is calculated assuming a constant heat flux $q_b$ through the bottom of the domain: $q_b=-6.2$ at ${{\textit{Re}}_{\tau ,0}}=170$ and $q_b=-19$ at ${{\textit{Re}}_{\tau ,0}}=636$ , in dimensionless units. These values were obtained directly from the DNS data. The tabulated values of $\mathcal{T}_{\textit {steady}}$ allow to estimate the computational cost needed to reach the equilibrium ice thickness. For this estimation, we took as reference a petascale Tier-0 supercomputer equipped with four GPUs per node and assumed perfect weak scalability (which is indeed what we observed within the range of MPI tasks employed for each simulation). Note that the computational cost is relative to the cheapest simulation (S1) and one computational cost unit is roughly equivalent to 50 GPU hours. We remark here that reaching the equilibrium ice thickness does not necessarily imply reaching a statistically steady state for the ice morphology, which is what we want to examine. For this reason, simulation runs were conducted beyond the equilibrium condition, implying that their effective cost is larger than that presented in table 1.

4.1. Morphodynamics in the subcritical regime

Figure 3 shows the ice layer for ${{\textit{Re}}_{\tau ,0}}=170$ and ${St}=1$ at statistically steady state. The layer is flipped upside down to better visualize the ice–water interface. In this case, the ice morphology is characterized by crests and valleys oriented in the streamwise direction (referred to hereinafter as turbulence-driven ice streaks) and separated by a characteristic wavelength $\lambda _y$ . Figure 4 shows the effect of the Stefan number on the ice morphology, more specifically on the fluctuating component of the ice thickness, $\xi ' (x,y,t)=\xi (x,y,t) -\langle \xi \rangle$ , is plotted for ${St}=0.1$ (figure 4 a), ${St}=1$ (figure 4 b) and ${St}=10$ (figure 4 c) at steady state. All the panels refer to the same rescaled time $t^*$ , defined in § 2. As it will be clarified later, the main features of the melting dynamics scale consistently with $t^*$ . The turbulence-driven ice streaks are clearly visible over the range of St we considered and their topological features do not change much with St. As reported by Perissutti et al. (Reference Perissutti, Marchioli and Soldati2024), these structures are directly generated by turbulence and should not in principle be affected by a change in the latent heat flux. The non-uniform heat flux provided by the near-wall fluid velocity streaks allows the ice crests and valleys to mark their footprint on the ice–water interface. However, some quantitative differences among the three cases can be observed. First, the maximum amplitude of the ice thickness fluctuations, $\xi '$ , decreases as St increases and the interface appears to be somewhat smoother in the streamwise direction, with $\xi '$ that remains almost unchanged along the $x$ direction. Conversely, at low St, the irregularity of the interface becomes stronger and changes in $\xi '$ become noticeable also along the streamwise direction. This behaviour is due to the fact that a larger latent heat contribution (i.e. larger St) reduces the volume of ice that is melted (or frozen) locally by the heat flux fluctuations induced by turbulence, reducing in turn the amplitude of $\xi '$ . This also has a time-filtering effect on the modifications of the ice thickness: a short pulse of high heat flux is less impactful on the final ice morphology and allows the latter to become more regular. The crest-to-trough height of the ice streaks ranges between roughly 5 and 15 w.u., a value large enough to allow dynamic interactions with the turbulent structures of the flow. These interactions cause the streamwise-oriented ice streaks to meander sideways, in a way that is similar to what the near-wall fluid velocity streaks do (Marchioli & Soldati Reference Marchioli and Soldati2002), albeit to a lesser extent. This can be explained considering that turbulent streaks are free to move sideways near a non-evolving flat wall, not equally so in the presence of ice streaks. Low-speed streaks of fluid tend to remain linked to the ice crests, while high-speed streaks tend to form more often in the troughs. This creates a self-reinforcing mechanism that pins the turbulent streaks by hindering their meandering and sustains the ice morphodynamics by increasing the clearance between crests and troughs. The correlation between ice crests and low-speed streaks (or equivalently, valleys and high-speed streaks) becomes stronger at lower St, because the typical time scales with which the interface evolves become shorter and hence closer to those of the fluid velocity streaks. The stronger the correlation, the larger the height difference between crests and troughs, explaining why they become more prominent at low St.

Figure 3.

Visualization of the ice layer, flipped upside down to highlight the interface morphology, in the subcritical regime ( ${{\textit{Re}}_{\tau ,0}} = 170$ ) and ${St} = 1$ . In this regime, the ice morphology is characterized by turbulence-driven ice streaks that are aligned along the streamwise direction and separated by a characteristic wavelength $\lambda _y$ .

Figure 4.

Influence of St on ice morphology in the subcritical regime ( ${{\textit{Re}}_{\tau ,0}} = 170$ ). The spatial distribution of the ice thickness fluctuation, $\xi '=\xi -\langle \xi \rangle$ , at steady state is shown for ${St}=0.1$ (a), ${St}=1$ (b) and ${St}=10$ (c). The snapshots refer to the same rescaled time ( $t^*=0.17$ ). Thicker ice regions are shown in white, thinner ice regions in blue. The ice morphology is not dramatically affected by St, yet some differences among the three cases can be noticed. At low St, the ice morphology is characterized by slight variations of $\xi '$ along the streamwise direction, these variations being damped at high St.

Figure 5 shows the evolution of the average ice thickness, $\langle \xi \rangle$ , as a function of the rescaled time, $t^*$ . The model discussed in § 2 predicts that, with such rescaling, the melting process should evolve independently of St and should follow the behaviour indicated by the black dashed line, which corresponds to the analytical solution of (2.13) that can be obtained when $q_b$ is provided. Here, we used $q_b = -6.2$ , a value obtained directly from the simulations. For ${St}=1$ and ${St}=10$ , the model captures well the evolution of the average ice thickness, while the prediction is just slightly worse for ${St}=0.1$ . Note that this value is rather close to the lower bound for the range of Stefan numbers within which the model is expected to perform reliably. The lower bound can be estimated using the condition (2.8), which yields ${St} \gg 0.05$ for the subcritical case, based on the initial value of the ice thickness. In spite of this, the model prediction is still fairly accurate, as also discussed in Appendix B. At the lower Stefan number, melting takes place very quickly in the early stages, while progressively slowing as the equilibrium condition, $\langle \xi \rangle =\langle \xi _{\textit{eq}}\rangle$ , is approached (at time $\mathcal{T}^*_{\textit {steady}} \approx 0.134$ ). This behaviour can be understood considering that the effects of heat diffusivity are more important at low St, and melting is mainly counteracted by the conductive heat flux through the ice layer. This heat flux is weak at the beginning of the simulation, when the ice layer is thick, but becomes stronger as ice gets thinner on a time scale that is dictated by the thermal diffusivity. At larger Stefan numbers ( ${St}=1$ and ${St}=10$ ), the latent heat flux is strong enough to overcome the effect of thermal diffusivity (which is neglected in the model) and produce a slower and more gradual melting, in better agreement with the theoretical prediction. The inset of figure 5 shows the evolution of $\langle \xi \rangle$ without rescaling: time is normalized by the eddy turnover time $\mathcal{T}_{\textit {eddy}}=h/u_\tau$ to highlight the dependence of the melting time scale on St. The simulation at ${St}=10$ is the one that takes longer to reach the steady state for $\langle \xi \rangle$ (roughly 300 eddy turnover times) and, hence, is the one that is computationally more expensive: the simulated time span covers over 600 eddy turnover times, almost three times that of the ${St}=1$ simulation and almost 10 times that of the ${St}=0.1$ simulation.

Figure 5.

Time evolution of the average ice thickness, $\langle \xi \rangle$ , at ${{\textit{Re}}_{\tau ,0}} = 170$ for ${St} = 0.1$ (orange line), ${St} = 1$ (blue line) and ${St} = 10$ (purple line) as a function of the rescaled time $t^*$ . The prediction from the analytical model is also shown (black dashed line). The model predicts a rescaled steady-state time equal to $\mathcal{T}^*_{\textit {steady}} \approx 0.134$ , also shown in the figure. The inset shows the evolution of the ice thickness when time is normalized by the eddy turnover time of the flow, $\mathcal{T}_{\textit {eddy}}$ , rather than rescaled using ${\textit{Pe}}$ and St.

Figure 6.

Time evolution of the ice fluctuation amplitude (i.e. the standard deviation of the ice thickness) $\langle \xi '^{2}\rangle ^{1/2}$ at ${{\textit{Re}}_{\tau ,0}} = 170$ for ${St} = 0.1$ (orange line), ${St} = 1$ (blue line) and ${St} = 10$ (purple line). In all cases, a statistically steady state is reached eventually, albeit at different times, ranging from 0.018 to 0.036 depending on St (as indicated by the light blue vertical band in the figure). This range defines the characteristic time scale of the turbulence-driven ice streaks formation, $\mathcal{T}^*_{\textit {streaks}}$ , in the subcritical regime. To allow a more direct comparison between the different time scales at play, in the following we take $\mathcal{T}^*_{\textit { streaks}}\approx 0.027$ , which is the average value among the different cases. Note that, even at steady state, $\langle \xi '^{2}\rangle ^{1/2}$ oscillates significantly around the averages steady-state value (marked for each case with a dashed line of the corresponding colour). For ${St} = 10$ , the standard deviation is significantly lower ( $\langle \xi '^{2}\rangle ^{1/2}\approx 0.0052$ ) compared with the other two cases ( $\langle \xi '^{2}\rangle ^{1/2}\approx 0.0128$ for ${St}=1$ and $\langle \xi '^{2}\rangle ^{1/2}\approx 0.0143$ for ${St}=0.1$ ).

Although the average ice thickness reaches an equilibrium condition, the local ice thickness continues to change over time while remaining in a statistically steady state. This is shown in figure 6, where the instantaneous standard deviation of the ice thickness, $\langle \xi '^{2}\rangle ^{1/2}$ , is plotted over the rescaled time $t^*$ . After a short initial transient, the standard deviation oscillates around a value that remains roughly constant over time. Because the turbulence-driven ice streaks are the only structures observed at low ${\textit{Re}}_{\tau ,0}$ , the extent of the transient may be interpreted as the characteristic time required to form the streaks, $\mathcal{T}^*_{\textit { streaks}}$ . Here, we assume that the transient ends as soon as the standard deviation first crosses its mean value at steady state. The crossover occurs at slightly different times, depending on St, within a range $0.018{-}0.036$ (with average value equal to 0.027) that is highlighted by the light blue vertical band in figure 6. This range defines the characteristic time scale of the turbulence-driven ice streaks formation, $\mathcal{T}^*_{\textit {streaks}}$ , in the subcritical regime, which is roughly five times shorter than the steady-state time $\mathcal{T}^*_{\textit {steady}}$ shown in figure 5. This indicates that the mechanism by which the ice streaks are generated, acts on time scales that are significantly shorter than those needed to reach the equilibrium condition. Figure 6 also confirms that the ice interface is less deformed at higher St, as observed qualitatively in figure 4. The average value of the standard deviation at ${St}=10$ is $\langle \xi '^{2}\rangle ^{1/2} \approx 0.0052$ , significantly lower than in the other two cases, which are fairly similar: $\langle \xi '^{2}\rangle ^{1/2} \approx 0.0128$ at ${St}=1$ and $\langle \xi '^{2}\rangle ^{1/2} \approx 0.0143$ at ${St}=0.1$ , respectively. This observation suggests that morphological changes of the interface are dominated by the latent heat effect at higher St, whereas the heat diffusivity is predominant at lower St. To further analyse the interface morphology, in figure 7 we show the spanwise spectra of the ice thickness, $|\hat {\xi }_y|$ . These spectra are obtained by averaging along the streamwise direction and over a short time interval $\Delta t^* = 0.034$ centred at time $t^*=0.17$ . They are plotted as a function of the spanwise wavenumber in w.u., $k_y^+ = k_y/{\textit{Re}}_{\tau }$ , with $k_y = 2\pi /\lambda _y$ . All spectra exhibit roughly the same decaying trend, with differences in amplitude that become smaller as St increases. This behaviour is consistent with previous observations and is confirmed by the inset of figure 7, which clearly shows that the spectra nicely overlap when rescaled by the Péclet and Stefan numbers as $|\hat {\xi }^*_y| ={\textit{Pe}} \sqrt {{St}+\mathcal{Q}} |\hat {\xi }_y|$ . This scaling relation was obtained empirically from the collapse of the rescaled spectra. Further investigation, for instance focused on the dynamical nature of the turbulence-driven ice streaks, is needed to provide a theoretical background for the ${\textit{Pe}} \sqrt {{St}+\mathcal{Q}}$ scaling.

Figure 7.

Spanwise spectra of the ice layer thickness, $|\hat {\xi }_y|$ , in the subcritical regime ( ${{\textit{Re}}_{\tau ,0}}=170$ ) for St = 0.1 (orange), ${St}=1$ (blue) and ${St}=10$ (purple). All the spectra are averaged along the streamwise direction and in time, over a short time interval $\Delta t^* = 0.034$ centred at time $t^*=0.17$ , taken at steady state. As shown in the inset, the spectra overlap almost perfectly when rescaled as $|\hat {\xi }^*_y| ={\textit{Pe}} \sqrt {{St}+\mathcal{Q}} |\hat {\xi }_y|$ .

4.2. Morphodynamics in the supercritical regime

The analysis carried out for the subcritical regime is applied here to the supercritical regime. In figure 8, we show an instantaneous visualization of the ice–water interface (flipped upside down to highlight its main morphological features) in the supercritical regime ( ${{\textit{Re}}_{\tau ,0}}=636$ ) at ${St}=1$ . First, it can be observed that the ice layer is thinner than in the subcritical regime, due to the highest heat flux resulting from the stronger turbulence of the water flow. In addition, the ice morphology appears more complex than in the subcritical regime due to the superposition of the turbulence-driven ice streaks with wave-like ripples of greater length scale that develop orthogonally to the flow direction. These structures are characterized by sharp crests separated by a wavelength, $\lambda _x$ , and exhibit a quasiperiodic behaviour along the streamwise direction. Ice ripples typically form as a result of a morphodynamic instability that emerges only at high Reynolds numbers. The sharpness of the crest, separated by concave depressions, is in line with what is observed in dissolution patterns and can be explained solely by geometric arguments stemming from the surface evolution (Chaigne et al. Reference Chaigne, Carpy, Massé, Derr, Courrech du Pont and Berhanu2023).

Figure 8.

Visualization of the ice layer (flipped upside down to highlight the interface morphology) in the supercritical regime ( ${{\textit{Re}}_{\tau ,0}} = 636$ ) at ${St} = 1$ . This regime is characterized by the presence of turbulence-driven ice streaks, having a characteristic wavelength $\lambda _y$ , and observed also in the subcritical regime, superposed to wavy ice ripples aligned with the spanwise flow direction, having a characteristic wavelength $\lambda _x$ .

Similar to the subcritical regime, the interface morphology is mildly affected by St, as shown in figure 9, where the fluctuating component of the ice thickness, $\xi '$ , is reported for ${St}=0.1$ (figure 9 a), ${St}=1$ (figure 9 b) and ${St}=10$ (figure 9 c). All plots refer to the same time instant, equal to the time at which the ice ripples are observed to reach their maximum height, $t^* \approx 0.0156$ . By visual inspection of figure 9, it is possible to appreciate that the main morphological features of the ice–water interface do not change dramatically with St: in all three cases, the ice–water interface is mainly characterized by the presence of ice ripples that have fairly similar wavelength and amplitude. The latter, in particular, can reach values of the order of $70$ – $80$ w.u., while the typical wavelength ranges between $\approx\! 530$ and $\approx\! 670$ w.u. These values are significantly lower than the minimum wavelength required for the onset of instability. This wavelength can be estimated using the models proposed in the literature, e.g. the model developed by Hanratty (Reference Hanratty1981), which predicts a minimum wavelength of approximately 1500 w.u., or the more recent model by Chedevergne et al. (Reference Chedevergne, Stuck, Olazabal-Loumé and Couzi2023), which predicts a minimum wavelength of approximately 1000 w.u. This discrepancy might be due to the finiteness of the computational domain in the wall-normal direction: the maximum size of the vortical structures is constrained by the channel height, $h$ , favouring the development of lower wavelengths. Nonlinear effects that the models are unable to capture may also explain the different values. Note that the streamwise length of the computational domain in the subcritical regime is equal to 1420 w.u. and is sufficient to accommodate at least two ripple wavelengths as measured in the supercritical regime. Considering that such wavelengths were found to be independent of the flow Reynolds number when expressed in wall units (Hsu et al. Reference Hsu, Locher and Kennedy1979; Thomas Reference Thomas1979), it seems reasonable to conclude that the lack of ice ripples in the subcritical regime is not related to the computational domain size.

Figure 9.

Influence of St on ice morphology in the supercritical regime ( ${{\textit{Re}}_{\tau ,0}} = 636$ ). The fluctuating component of the ice thickness map, $\xi '=\xi -\langle \xi \rangle$ , is shown for ${St}=0.1$ (a), ${St}=1$ (b) and ${St}=10$ (c). All plots are taken at the rescaled time $t^* \approx 0.0156$ , which is when the amplitude of the ice ripples is observed to become the highest. Thicker ice regions are shown in white, thinner ice regions in blue. For all St, the supercritical morphology is characterized by the presence of ice ripples, which are fairly similar in shape. The main effect of St is to reduce the size of the superposed turbulence-driven ice streaks, which are less visible at high St.

Figure 10.

Time evolution of the average ice thickness, $\langle \xi \rangle$ , at $ {{\textit{Re}}_{\tau ,0}} = 636$ for ${St}=0.1$ (orange), ${St}=1$ (blue) and ${St} = 10$ (purple). The prediction of the analytical model is represented by the black dashed line. The steady-state time predicted by the model is also shown: $\mathcal{T}^*_{\textit {steady}} \approx 0.023$ . Similarly to the subcritical case, the prediction is fairly accurate for large St, while becoming less precise for low St. The inset shows the evolution of the ice thickness when time is normalized by the eddy turnover time, $\mathcal{T}_{\textit {eddy}}$ , rather than rescaled using ${\textit{Pe}}$ and St.

Despite the limited St effects, figure 9 shows some interesting modifications: at high St, the ice ripples seem less blurred, most likely because the latent heat flux acts as a temporal filter and mitigates the effect of the heat fluctuations on the interface. Another reason may be that, while the amplitude of ice ripples seems unaffected by St, the depth of the overlapping turbulence-driven ice streaks does depend on St, as discussed in the previous section. This results in a smaller height of the turbulence-driven ice streaks relative to the size of the ice ripples, which appear more defined.

Figure 10 shows the evolution of the average ice layer thickness, $\langle \xi \rangle$ , as a function of the rescaled time, $t^*$ . The model prediction, obtained solving (2.13) with $q_b = -19$ as provided by the numerical simulation, is again indicated by the black dashed line. Similarly to the subcritical case, the accuracy of the model improves with increasing St: indeed, the simulation at ${St}=10$ (purple line) matches very well with the model prediction. The prediction for ${St}=1$ (blue line) is still fairly accurate, indicating that the effect of the latent heat becomes predominant: melting slows down and the evolution of $\langle \xi \rangle$ becomes linear in time. Differences become more significant at ${St}=0.1$ since, in this case, the condition (2.8) predicts that the model is accurate if ${St} \gg 1.35$ . This is evident especially during the early stages of the interface evolution, when fast melting occurs due to heat diffusivity effects and the error of the model is large because the interface is still far from equilibrium (see Appendix B). Similarly to the subcritical case, as the interface approaches equilibrium, the error decreases, allowing the model to provide reasonably accurate predictions even for ${St} = 0.1$ and ${St} = 1$ . Note that, for the case ${St}=1$ , the melting rate (given by the slope of the curves in figure 10) is larger than the value predicted by the model within the time interval from $t^* \simeq 0.01$ and $t^* \simeq 0.015$ . As it will be shown later, this is the time window within which the ice ripples reach their maximum amplitude. Interestingly, the model is always able to capture correctly the equilibrium ice thickness, $\langle \xi _{\textit{eq}}\rangle$ , which is independent of St and has the same value in all the simulations.

The rescaled steady-state time required to reach the equilibrium is also similar in all the three cases: $\mathcal{T}^*_{\textit {steady}}\approx 0.023$ , computed from (2.14). This implies that the model can capture with acceptable accuracy the relevant time scales of the melting process even in the supercritical regime. Note that, if no time rescaling is considered, then these time scales are much bigger for larger St, as shown in the inset of figure 10, where $\langle \xi \rangle$ is plotted as a function of the non-rescaled time $t$ normalized by the eddy turnover time, $\mathcal{T}_{\textit {eddy}}$ . The inset shows that, as in the subcritical regime (see discussion of figure 5), the most expensive simulation from a computational point of view is the one at ${St}=10$ , which covers almost 250 eddy turnover times.

A notable difference between the subcritical and supercritical regimes is that, in the latter regime, the standard deviation, $\langle \xi '^{2}\rangle ^{1/2}$ , does not simply oscillate around a constant, steady-state value but rather exhibits an initial growth followed by a decrease before reaching the steady state, as shown in figure 11. The fluctuation amplitude grows until a maximum value (highlighted by the dot in each curves) is reached. The maximum is reached at different times depending on St, but within a relatively narrow time range (marked with the light blue region in the figure), centred at $t^*\approx 0.0156$ . This value represents the characteristic time scale required to fully form the ripples (see also discussion of figure 9). Once the maximum amplitude is attained, $\langle \xi '^{2}\rangle ^{1/2}$ decreases until the steady state is reached (again, at different times depending on St). This behaviour can be understood considering that the fluctuation amplitude is influenced not only by the presence of the turbulence-driven ice streaks, as is the case in the subcritical regime, but also (mostly, in fact) by the ice ripples, which are much larger in size and hence more impactful on the behaviour of $\langle \xi '^{2}\rangle ^{1/2}$ . As will be detailed in § 4.3, we can anticipate that ice ripples undergo the same dynamics (initial growth, maximum amplitude, subsequent decay), suggesting that the peak of $\langle \xi '^{2}\rangle ^{1/2}$ and the maximum amplitude of the ice ripples occur approximately at the same time, namely $\mathcal{T}^*_{\textit {ripples}}$ .

Figure 11.

Time evolution of the ice fluctuation amplitude (equal to the standard deviation of the ice thickness), $\langle \xi '^{2}\rangle ^{1/2}$ , in the supercritical regime ( ${{\textit{Re}}_{\tau ,0}} = 636$ ) for ${St} = 0.1$ (orange), ${St} = 1$ (blue) and ${St} = 10$ (purple). Due to the presence of the ice ripples, $\langle \xi '^{2}\rangle ^{1/2}$ grows until reaching a maximum value (marked by a dot for each case). The maximum value of $\langle \xi '^{2}\rangle ^{1/2}$ is attained within a time window (highlighted by the light blue vertical band) centred at $t^*\approx 0.0156$ . This value represents the characteristic time scale required to fully form the ripples and is indicated as $\mathcal{T}^*_{\textit {ripples}}$ in the figure. Once the maximum amplitude is attained, $\langle \xi '^{2}\rangle ^{1/2}$ decreases until the steady state is reached (at different times depending on St). The time required to form the turbulence-driven ice streaks (also shown) is much shorter than the ice ripples time scale: $\mathcal{T}^*_{\textit {streaks}}\approx 8.14 \boldsymbol{\cdot }10^{-4}$ .

The trends shown in figure 11 do not allow a straightforward quantification of the time scale required to form the turbulence-driven ice streaks, $\mathcal{T}^*_{\textit {streaks}}$ , in the supercritical regime. Due to the presence of the ripples, $\mathcal{T}^*_{\textit {streaks}}$ cannot be simply set equal to the steady-state time for $\langle \xi '^{2}\rangle ^{1/2}$ , as done for the subcritical regime, see figure 6. The reason is that the steady state is reached only when the ice ripples stabilize. The alternative definition we adopted in this case is to set $\mathcal{T}^*_{\textit {streaks}}$ equal to the steady-state time for the spanwise-spectrum amplitude of the ice interface, $|\hat {\xi }_y|$ . This spectrum is evaluated at the wavenumber that corresponds to the typical wavelength of the turbulence-driven ice streaks ( $\lambda _y^+\approx 100$ ), not shown here for brevity. More specifically, for each case, $\mathcal{T}^*_{\textit {streaks}}$ , is set equal to the rescaled time at which $|\hat {\xi }_y|_{\lambda ^+=100}$ first falls within 5 % of its time-averaged value (computed from $\mathcal{T}^*_{\textit {streaks}}$ onward). This results in a narrow range of values of $\mathcal{T}^*_{\textit {streaks}}$ , centred around $8.14\times 10^{-4}$ , between the different simulations, highlighted by the light blue vertical band in figure 11. In the supercritical case, the difference between $\mathcal{T}^*_{\textit {steady}}$ and $\mathcal{T}^*_{\textit {streaks}}$ , corresponding to a ratio $\mathcal{T}^*_{\textit {steady}} / \mathcal{T}^*_{\textit {streaks}} \approx 30$ , is larger than in the subcritical case, where corresponding to a ratio $\mathcal{T}^*_{\textit {steady}} / \mathcal{T}^*_{\textit {streaks}} \approx 5$ . This is likely due to the size of the turbulence-driven ice streaks, which is significantly smaller than the domain size (which increases linearly with ${\textit{Re}}_{\tau ,0}$ in each flow direction) in the supercritical case. The time scale of the turbulence-driven ice streaks is also much smaller than that of the ice ripples, $\mathcal{T}^*_{\textit {ripples}}$ , their ratio being $\mathcal{T}^*_{\textit {ripples}} / \mathcal{T}^*_{\textit {streaks}} \approx 20$ . This indicates that the two phenomena are well separated in time and are governed by independent physical mechanisms. We also notice that the steady-state value of $\langle \xi '^{2}\rangle ^{1/2}$ is attained between $t^* \approx 0.023$ and $t^* \approx 0.025$ , which is approximately the same time at which the steady state for $\langle \xi \rangle$ is reached (see value of $\mathcal{T}^*_{\textit {steady}}$ in figure 10).

As mentioned, it is not possible to filter out the effect of the ice ripples on the ice thickness fluctuations from figure 11. To isolate and characterize the effect of the turbulence-driven ice streaks, we thus examine the spanwise spectra shown in figure 13. These spectra are unaffected by the presence of the ripples, which induce significant interface deformations along the streamwise direction only. As done in figure 7, spectra are averaged in space (along $x$ ) and time, over a time interval $\Delta t^*=0.003$ centred at $t^*=0.0156 \equiv \mathcal{T}^*_{\textit {ripples}}$ : within this interval, the ice ripples reach their maximum amplitude regardless of St. Figure 13 shows clearly that, also in the supercritical regime, spanwise spectra share the same decaying trend and overlap nicely when rescaled as $|\hat {\xi }_y^*| = |\hat {\xi }_y|{\textit{Pe}} \sqrt {{St} + \mathcal{Q}}$ (see figure inset). Additionally, the inset shows the rescaled spectra for the subcritical cases (dashed lines), which collapse on top of the simulations at ${{\textit{Re}}_{\tau ,0}}=636$ . This indicates that the morphology of the turbulence-driven ice streaks is similar over the range of St investigated in our study, but also over the range of ${\textit{Re}}_{\tau ,0}$ (namely ${\textit{Pe}}$ ) characterizing the different regimes. The emergence of ripples does not affect the spanwise features of the ice–water interface, but increases the variability of the interface morphology along the streamwise direction, as confirmed by the streamwise spectra, $|\hat {\xi }_x|$ , shown in figure 12. These spectra are averaged in space, along the spanwise direction, and in time, over the same interval of figure 13. Three peaks, indicating the presence of the ice ripples, are clearly visible at slightly different values of the dimensionless streamwise wavenumber, $k_x$ : more precisely, $k_{x,max}=8$ , $k_{x,max}=10$ and $k_{x,max}=9$ for ${St}=0.1$ , ${St}=1$ and ${St}=10$ , respectively. Note that no peak is observed in the subcritical spectra ( ${{\textit{Re}}_{\tau ,0}} = 170$ ), shown with dashed lines in the inset of figure 12. Recalling that $k_x^+ = k_x/{\textit{Re}}_{\tau }=2\pi /\lambda _x^+$ and $L_x^+ = 2 \pi {\textit{Re}}_{\tau }$ , it follows that the supercritical peak values of $|\hat {\xi }_x|$ correspond to characteristic wavelengths of the ripples $\lambda _x^+ = L_x^+/8$ , $L_x^+/10$ and $L_x^+/9$ . Once rescaled, these wavenumbers cover a narrow range of values centred at $k_x^+ = k_x/{\textit{Re}}_{\tau } \approx 10^{-2}$ : this implies that $\lambda _x$ is weakly affected by St, confirming the qualitative observations drawn from figure 9. Also poorly affected by the Stefan number is the amplitude of the ice ripples, given that the peak values of $|\hat {\xi }_y|$ are almost identical in the three simulations: $|\hat {\xi }_x|_{max} \approx 3 \times 10^{-3}$ .

Figure 12.

Streamwise spectra of the ice thickness, $|\hat {\xi }_x|$ , in the supercritical regime ( ${{\textit{Re}}_{\tau ,0}} = 636$ , solid lines) for ${St} = 0.1$ (orange), ${St} = 1$ (blue) and ${St} = 10$ (purple). The spectra are averaged both in space, along the spanwise direction, and time, within the interval $\Delta t^*=0.003$ centred at time $t^*=0.0156 \equiv \mathcal{T}^*_{\textit {ripples}}$ . In all cases, spectra are characterized by a peak, indicated with a filled circle, that occurs at slightly different wavenumbers $k_{x,max}$ , and is associated with the presence of the ice ripples. Spectra in the inset, corresponding to the ${{\textit{Re}}_{\tau ,0}} = 170$ case, exhibit no peak.

Figure 13.

Spanwise spectra of the ice layer thickness, $|\hat {\xi }_y|$ , in the supercritical regime ( ${{\textit{Re}}_{\tau ,0}}=636$ ) for St = 0.1 (orange), ${St}=1$ (blue) and ${St}=10$ (purple). All the spectra are averaged along the streamwise direction and in time, over a short time interval $\Delta t^* = 0.003$ centred at time $t^*=0.0156 \equiv \mathcal{T}^*_{\textit {ripples}}$ . The inset shows the rescaled spectra, $|\hat {\xi }^*_y| ={\textit{Pe}} \sqrt {{St}+\mathcal{Q}} |\hat {\xi }_y|$ , also for the subcritical regime (dashed lines). All curves overlap, indicating that the ( ${\textit{Pe}}$ , St) scaling holds regardless of the morphodynamic regime.

4.3. Characterization of ripple dynamics

As discussed in the previous section, the key morphological feature in the supercritical regime is the presence of the ice ripples. In this section, we examine their growth over time, which occurs because the heat flux supplied by the warm water stream is not in phase with the ice layer thickness, $\xi$ . In particular, ice ripples can grow only if the heat flux lags behind $\xi$ by an angle greater than $\pi /2$ (Ashton & Kennedy Reference Ashton and Kennedy1972; Gilpin et al. Reference Gilpin, Hirata and Cheng1980). Physically, this means that the amplitude of the ripples increases only if the heat flux perturbations induced by the wavy interface reach a minimum near the crest and a maximum near the trough.

To quantify the growth of the ice ripples and its dependence on St, in figure 14 we show the time evolution of $|\hat {\xi }_x|_{max}$ , which corresponds to the peak value of the spectra in figure 12 and represents the amplitude of the frequency mode associated with the ice ripples. It is apparent that the qualitative behaviour of $|\hat {\xi }_x|_{max}$ does not depend on St. The amplitude of the ripples grows in the early stages of the simulations, during the transient dominated by the already-formed turbulence-driven ice streaks. Around time $t^* \approx 0.0085$ , namely at approximately one-third of the supercritical, steady-state time $\mathcal{T}^*_{\textit {steady}}$ , the mode associated with the ice ripples becomes dominant in the streamwise spectra, thus marking the transition to a phase in which the ice morphology is dominated by the ice ripples growth. The ice ripples reach their maximum amplitude around $t^* \approx 0.0160$ for ${St}=0.1$ , $t^* \approx 0.0149$ for ${St}=1$ and $t^* \approx 0.0162$ for ${St}=10$ , namely within a relatively narrow interval (highlighted in light blue in figure 14) centred at $\mathcal{T}^*_{\textit {ripples}} \approx 0.0156$ . This confirms the large time separation that characterizes the formation of streaks and ripples ( $\mathcal{T}^*_{\textit {ripples}} / \mathcal{T}^*_{\textit {streaks}} \gg 1$ ), and suggests that these two processes can be decoupled and modelled separately.

Figure 14.

Time evolution of the maximum value, $|\hat {\xi }_x|_{max}$ , of the streamwise spectrum of the ice thickness as a function of the rescaled time, $t^*$ , at varying St. As shown in figure 12, the peak value is obtained at $k_{x} = 8$ for ${St}=0.1$ (orange), at $k_{x}=10$ for ${St}=1$ (blue) and at $k_x=9$ for ${St}=10$ (purple). Regardless of St, the amplitude grows in the early stages of the simulation, reaching a maximum (marked by a dot) within a relatively narrow range, represented by the larger light blue vertical band. This range does depend on St and is centred at $t^* \approx 0.0156 \equiv \mathcal{T}^*_{\textit {ripples}}$ . Note that $\mathcal{T}^*_{\textit {ripples}}$ is significantly larger than $\mathcal{T}^*_{\textit {streaks}}$ while being of the same order of $\mathcal{T}^*_{\textit {steady}}$ . The final stages of the evolution are characterized by a decay of the amplitude due to the increasing strength of heat conduction within the ice layer, which eventually leads to the steady state.

We observe also that $\mathcal{T}^*_{\textit {ripples}}$ is roughly equal to $2/3$ of the steady-state time reported in figure 10, namely $\mathcal{T}^*_{\textit {ripples}} / \mathcal{T}^*_{\textit {steady}} \simeq {O}(1)$ : this implies that the morphodynamic evolution of the ice ripples is strongly dependent on the average thickness of the ice layer. This is consistent with the theory formulated by Ashton & Kennedy (Reference Ashton and Kennedy1972). According to this theory, the growth rate of the ice ripples depends on two competing effects: the heat flux supplied by the warm water flow, which increases the amplitude of the ripples, and the diffusive heat flux from the ice layer, which flattens the ripples. The latter is closely related to the average thickness of the ice layer since, as $\langle \xi \rangle$ decreases, the diffusive heat flux through the ice layer increases (in the quasistationary approximation, it is inversely proportional to $\langle \xi \rangle$ ). At some point, the diffusive heat flux overcomes the effect of the heat flux supplied by the water stream, leading to a reduction of $|\hat {\xi }_x|_{max}$ and explaining the decaying trend observed in figure 10 at times larger than $\mathcal{T}^*_{\textit {ripples}}$ . The decay suggests that ripples would not form if the supercritical simulations were started using the ice layer morphology at equilibrium as initial condition. This shows that the formation of ice ripples is favoured when the ice layer is shrinking. The maximum value of the ice ripples amplitude for ${St}=1$ and ${St}=10$ is similar, while being larger for ${St}=0.1$ . This could be due to the fact that the diffusive heat flux through the ice layer does not immediately adjust to changes in the ice morphology but is delayed due to the time scales at which the thermal diffusivity acts. This delay may hamper the capability of heat diffusion within the ice to flatten the interface. Besides this low-St effect, however, the limited variability we observe for the maximum amplitudes is not surprising. It is reasonable to assume that heat flux anomalies at the ice–water interface are mostly influenced by the ice morphology and therefore do not depend on the latent heat.

Figure 14 also shows a nearly exponential growth of $|\hat {\xi }_x|_{max}$ for ${St}=1$ and ${St}=10$ . This is again consistent with the theory by Ashton & Kennedy (Reference Ashton and Kennedy1972), which associates the growth of amplitude of the ripples to a time scale that is directly proportional to the thermal diffusivity (namely ${\textit{Pe}}^{-1}$ ) and inversely proportional to the latent heat (namely ${St}^{-1}$ ). As a result, the growth rate of the ripples at larger St should be similar when evaluated in terms of the rescaled time $t^*$ . We remark here that, in our problem, $t^*$ also depends on the additional sensible heat contribution $\mathcal{Q}$ , defined in § 2.2. This contribution is not accounted for in the theory by Ashton & Kennedy (Reference Ashton and Kennedy1972), which is based on the quasistationary approximation for the temperature fields. As discussed in § 2, this condition holds only for very large St, so the effect of $\mathcal{Q}$ must be included for a proper comparison at low St, when such approximation does not hold. For ${St}=0.1$ , we observe an earlier onset of the growth due to the low latent heat flux, quickly followed by a slower increase compared with the two cases at higher St, limited by the propagation speed of thermal diffusivity.