2. Underlying equations, and the control and response parameters

In the present study, we assume the flow under consideration to be incompressible, such that the density field $\rho$ is determined by the temperature ( $T ( \boldsymbol{x},t )$ ) and salinity ( $S ( \boldsymbol{x},t )$ ) fields. The governing equations for our direct numerical simulation, incorporating the phase-field variable $\phi$ needed for melting, along with the incompressibility condition $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0$ , read as follows:

(2.1) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} &= -\frac {\boldsymbol{\nabla }p}{\rho _{0}}+\nu \boldsymbol{\nabla }^{2}\boldsymbol{u}-\frac {g{\rho }'}{\rho _{0}} \boldsymbol{e}_{y}-\frac {\phi \boldsymbol{u}}{\eta }, \end{align}

(2.2) \begin{align} \frac {\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }T &= \kappa _{T}\boldsymbol{\nabla }^{2}T+\frac {\mathcal{L}}{c_{\!p}}\frac {\partial \phi }{\partial t}, \end{align}

(2.3) \begin{align} \frac {\partial S}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }S &= \kappa _{S}\boldsymbol{\nabla }^{2}S+\frac {1}{1-\phi +\delta }\left ( S\frac {\partial \phi }{\partial t}- \kappa _{S}\boldsymbol{\nabla }\phi \boldsymbol{\cdot }\boldsymbol{\nabla }S\right)\!, \end{align}

(2.4) \begin{align} \frac {\partial \phi }{\partial t} &= C\boldsymbol{\nabla }^{2}\phi -\frac {C}{\varepsilon ^{2}}\phi \left ( 1-\phi \right )\left ( 1-2\phi +\frac {\varepsilon }{\gamma }\left ( T+mS \right )\right)\!. \end{align}

Here, $\boldsymbol{u} ( \boldsymbol{x},t )$ is the velocity field, ${\rho }'=\rho -\rho _{0}$ is the fluctuating density from a reference value $\rho _{0}$ , $p$ is the kinematic pressure, $\nu$ is the kinematic viscosity of the liquid, $g$ is the gravitational acceleration in the vertical direction $\boldsymbol{e}_{y}$ , and $\kappa _{T}$ and $\kappa _{S}$ are the diffusivities of temperature and salinity, respectively. In simulations, we employ the phase-field method (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020) to simulate the melting process of a solid in a multi-component fluid, which has been extensively used in numerous prior numerical investigations of melting ice (Favier et al. Reference Favier, Purseed and Duchemin2019; Couston et al. Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2021; Hester et al. Reference Hester, McConnochie, Cenedese, Couston and Vasil2021; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023a , Reference Yang, Howland, Liu, Verzicco and Lohseb ). In this method, the phase field variable $\phi$ is integrated in time and space, and smoothly transitions from a value of 1 in the solid to a value of 0 in the liquid, with the interface located at $\phi =1/2$ . Here, $\mathcal{L}$ is the latent heat of fusion and $c_{\!p}$ is the specific heat capacity. The liquidus slope $m={0.056}\,\textrm {K} ({\textrm {g}\,\textrm {kg}^{-1}})^{-1}$ reflects how the presence of salinity at the ice–water interface lowers the local melting temperature.

The phase-field model employed in this study follows the second-order formulation of Hester et al. (Reference Hester, Couston, Favier, Burns and Vasil2020) and properly reflects the Gibbs–Thomson effect (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020). The parameters $\varepsilon$ , $\eta$ , $C$ , $\gamma$ and $\delta$ all refer to the phase-field model, which can be explained as follows. The parameter $\varepsilon$ is used to measure the diffuse interface thickness, which is set to be equal to the grid spacing following the convergence test of Favier et al. (Reference Favier, Purseed and Duchemin2019). The limit $\varepsilon \rightarrow 0$ leads to the exact Stefan boundary conditions for $S$ and $T$ at the liquid–solid interface:

(2.5) \begin{align} S^{\left ( l \right )}u_{n} &= -\, \kappa _{S} \frac {\partial S^{\left ( l \right )}}{\partial n}, \end{align}

(2.6) \begin{align} \mathcal{L}u_{n} &= c_{\!p}\kappa _{T}\left ( \frac {\partial T^{\left ( s \right )}}{\partial n} - \frac {\partial T^{\left ( l \right )}}{\partial n} \right)\!, \end{align}

where $u_{n}$ is the normal velocity of the interface between the solid and the liquid phases, $n$ represents the normal direction of the interface, and the superscripts $ (s )$ and $ (l )$ represent the solid and liquid phases, respectively. The penalty parameter $\eta$ is used to decay the velocity to zero in the solid phase and its value is set to be equal to the time interval (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020). Furthermore, a direct forcing method is applied to set the velocity to zero for $\phi \gt 0.9$ to avoid spurious motions in the solid phase (Howland Reference Howland2022). The diffusivity $C$ of the phase field equation (2.4) is defined by $C=6\gamma \kappa _{T}/ (5\varepsilon \mathcal{L} )$ , where $\gamma$ is the surface energy coefficient related to the Gibbs–Thompson effect (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020; Howland Reference Howland2022; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023a ). We set $C=1.2\kappa _T$ and $\varepsilon \Delta T/\gamma =10$ , where $\Delta T$ is the initial temperature difference between the solid and liquid phases. The small parameter $\delta \ll 1$ is solely used to stabilise the terms on the right-hand side of the salinity equation (2.3) (Howland Reference Howland2022). Further description and validation of the phase-field method can be found from Hester et al. (Reference Hester, Couston, Favier, Burns and Vasil2020) and Howland (Reference Howland2022).

In the simulations, the initial temperature of the water surrounding the ice is set to $T_{\infty }={20}\,^\circ {\textrm {C}}$ . Since this is much larger than the maximum-density temperature $T_c={4}\,^\circ {\textrm {C}}$ , the Oberbeck–Boussinesq approximation is employed, ignoring the density anomaly that occurs at ${4}\,^\circ {\textrm {C}}$ for water. The fluctuating density is assumed to depend linearly on two scalar fields, namely the temperature $T$ and salinity $S$ ,

(2.7) \begin{align} \rho '=\rho _0 \left ( \beta _{S}S - \beta _{T}T \right ), \end{align}

where $\beta _T$ is the thermal expansion coefficient and $\beta _S$ is the haline contraction coefficient. We verified that the results of simulations considering the density anomaly are nearly identical to those of simulations ignoring the density anomaly. Obviously, this does no longer hold for an ambient water temperature $T_{\infty } \approx {4}\,^\circ {\textrm {C}}$ .

The system can be controlled using five dimensionless control parameters, namely the Rayleigh numbers of temperature and salinity, the Prandtl number, the Schmidt number and the Stefan number:

(2.8) \begin{align} \textit{Ra}_{T}=\frac {g\beta _{T} \Delta TH^{3}}{\nu \kappa _{T}},\quad \textit{Ra}_{S}=\frac {g\beta _{S}\Delta SH^{3}}{\nu \kappa _{S}},\quad \textit{Pr}=\frac {\nu }{\kappa _{T}},\quad {Sc}=\frac {\nu }{\kappa _{S}},\quad {Ste}=\frac {c_{\!p}\Delta T }{\mathcal{L}}, \end{align}

where $\Delta T=T_{\infty }-T_{s}$ and $\Delta S=S_{\infty }-S_{s}$ are the temperature and salinity differences between the surrounding water and the ice, respectively. Here, $T_s={0}\,^\circ {\textrm {C}}$ is the temperature of the ice, which is the equilibrium melting temperature for freshwater. Furthermore, $S_{\infty }$ and $S_s={0}\,{\textrm {g}\,\textrm {kg}^{-1}}$ are the salinity of the ambient water and the ice, respectively.

The Lewis number ${Le}$ and the density ratio $R_\rho$ corresponding to (2.8) are given by

(2.9) \begin{equation} {Le}=\frac {\kappa _{T}}{\kappa _{S}}=\frac {{Sc}}{\textit{Pr}}, \quad R_\rho \equiv \frac {\textit{Ra}_S \textit{Pr}}{\textit{Ra}_T {Sc}} = \frac {\beta _S \Delta S}{\beta _T \Delta T}, \end{equation}

Here, the density ratio $R_\rho$ is introduced to quantify the ratio between the strengths of the saline and thermal buoyancies.

While the experiments and simulations are initialised using the control parameters $\textit{Ra}_T$ and $\textit{Ra}_S$ , the system is physically driven by the differences in temperature and salinity between the ambient water and the ice–water interface, which are $\Delta T_d = T_\infty - T_i$ and $\Delta S_d = S_\infty - S_i$ , respectively. Here, $T_{i}$ and $S_{i}$ denote the temperature and salinity at the ice–water interface. Accordingly, the corresponding Rayleigh numbers are given by

(2.10) \begin{align} \textit{Ra}_{T_d}=\frac {g\beta _{T} \Delta T_{d}H^{3}}{\nu \kappa _{T}},\quad \textit{Ra}_{S_d}=\frac {g\beta _{S}\Delta S_{d}H^{3}}{\nu \kappa _{S}}. \end{align}

Since $T_i$ and $S_i$ are actively changing during the experiments and simulations, $\textit{Ra}_{T_d}$ and $\textit{Ra}_{S_d}$ are response parameters, as opposed to the control parameters $\textit{Ra}_T$ and $\textit{Ra}_{S}$ .

For the experiments, both $T_i$ and $S_i$ are assumed to be zero, which is a reasonable assumption for the relatively high ambient temperature considered in this study ( $T_\infty \approx {20}\,^\circ {\textrm {C}}$ ). In this case, Josberger & Martin (Reference Josberger and Martin1981) experimentally verified that the interface salinity is less than 1 % of the ambient salinity at least up to $S_\infty =$ ${35}\,{\textrm {g}\,\textrm {kg}^{-1}}$ . According to the liquidus condition, the interface temperature is therefore less than 0.1 % of the ambient temperature. Note that the interface salinity becomes more relevant as the ambient temperature approaches the freezing temperature (Kerr & McConnochie Reference Kerr and McConnochie2015).

For the simulations, the interface temperature $T_i$ and salinity $S_i$ are computed as time-averaged surface-integral values over the ice–water interface. Specifically, $S_i$ is computed by $ S_i= \langle \int \left |\boldsymbol{\nabla }\phi \right |S\,{\textrm d}\varOmega /\int \left |\boldsymbol{\nabla }\phi \right |\,{\textrm d}\varOmega \rangle _t$ , where $\varOmega$ represents the computational domain. The interface temperature $T_i$ is subsequently determined from $S_i$ via the liquidus condition $T_i=T_{m}-mS_{i}$ , where $T_{m}={0}\,^\circ {\textrm {C}}$ is the melting point of fresh water.

In this paper, the normalised mean melt rate $\widetilde {f}$ is defined by $\widetilde {f}=t_{D}/t_{m}$ , where $t_{m}$ represents the time needed to melt $V_{m}=70\,\%$ of the initial volume in three-dimensional (3-D) simulations and experiments and of the initial area in two-dimensional (2-D) simulations. Here, $t_{D}=H^{2}/\kappa _{T}$ is the diffusion time scale of temperature. We also tried to use different values of $V_{m}$ to calculate the melt time $t_{m}$ and the normalised mean melt rate $\widetilde {f}$ in 2-D simulations. This only resulted in a change of the absolute values of $t_{m}$ and $\widetilde {f}$ , but the trend and scaling remained the same. This invariance is further corroborated by the approximately constant instantaneous melt rate observed in the 2-D simulations, as discussed in Appendix C. Therefore, the value of $V_{m}$ has no impact on the qualitative results of the present study.

According to the Stefan boundary condition for salinity equation (2.5), the dimensionless instantaneous melt speed $\widetilde {u}_{n}$ can be related to the dimensionless instantaneous salt flux ${\textit{Nu}}_{S}$ :

(2.11) \begin{equation} \widetilde {u}_{n}=\frac {u_{n}}{U_{0}}=-\frac {1}{U_{0}}\frac {\kappa _{S}}{S}\frac {\partial S}{\partial n}=\frac {{\textit{Nu}}_{S}}{{Le}}\frac {S_{\infty }-S_{i}}{S_{i}}\frac {H}{\delta _S(t)}, \end{equation}

where $U_{0}=H/t_{D}$ is the thermal diffusion velocity scale, $\delta _S ( t )$ is the saline boundary layer thickness at time $t$ and ${\textit{Nu}}_{S}=- ( \partial S/\partial n )/ ( \Delta S_{d}/\delta _S(t) )$ . Since the mean melt rate $\widetilde {f}$ is per definition equal to the time-average of $\widetilde {u}_{n}$ , (2.11) expresses that $\widetilde {f}$ must be proportional to the dimensionless time-averaged salt flux $\overline {\textit{Nu}}_{S}$ .

Similar to (2.11), we relate the dimensionless instantaneous melt speed $\widetilde {u}_{n}$ to the dimensionless instantaneous heat flux $\textit{Nu}_{T}$ , according to the Stefan boundary condition for temperature equation (2.6):

(2.12) \begin{align} \widetilde {u}_{n}=\frac {u_{n}}{U_{0}}=-\frac {1}{U_{0}}\frac {\kappa _{T} C_{\!p}}{\mathcal{L}}\frac {\partial T}{\partial n}={Ste}\,\textit{Nu}_{T} \frac {T_{\infty }-T_{i}}{T_{\infty }} \frac {H}{\delta _T(t)}, \end{align}

where ${\textit{Nu}}_{T}=- ( \partial T/\partial n )/ ( \Delta T_{d}/\delta _T(t) )$ and $\delta _T ( t )$ is the thermal boundary layer thickness at time $t$ . Throughout the entire melting process, the ice temperature $T_{\textit{ice}}$ remains within the range of $-0.4\,^{\circ }\textrm{C}$ to $0\,^{\circ }\textrm{C}$ , indicating that the heat flux from the interior of the ice to the ice–water interface is much smaller than the heat flux from the surrounding water to the interface. Therefore, we assume that the term $\partial T^{ ( s )}/\partial n$ is much smaller compared with the term $\partial T^{ ( l )}/\partial n$ . Accordingly, the mean melt rate $\widetilde {f}$ should be proportional to the dimensionless time-averaged heat flux $\overline {\textit{Nu}}_{T}$ . Several previous studies (Josberger & Martin Reference Josberger and Martin1981; Kerr Reference Kerr1994; Kerr & McConnochie Reference Kerr and McConnochie2015; Howland, Verzicco & Lohse Reference Howland, Verzicco and Lohse2023) have reported that the ratio of thermal to solutal boundary layer thicknesses scales as $\delta _T(t)/\delta _S(t) \propto {Le}^{1/2}$ for Lewis numbers of order ${Le} = \mathcal{O}(100)$ . Adopting this scaling, we assume $\delta _T(t)/\delta _S(t) = {Le}^{1/2}$ , and we find the relation between $\overline {\textit{Nu}}_{S}$ and $\overline {\textit{Nu}}_{T}$ , namely (Holland & Jenkins Reference Holland and Jenkins1999; Couston Reference Couston2024)

(2.13) \begin{equation} \overline {\textit{Nu}}_{S}={Ste}\,{Le}^{1/2}\,\frac {S_{i}}{S_{\infty }-S_{i}}\frac {T_{\infty }-T_{i}}{T_{\infty }}\,\overline {\textit{Nu}}_{T}. \end{equation}

Figure 1.

(a) Schematic of the experimental set-up. A glass tank is filled with water. The cylinder is placed at the centre of the tank and is kept in place using a holder made from POM which is thermally poorly conducting. The ice has a hemispherical cap at the bottom, and is flat at the top. (b) and (c) Schematic of the numerical set-up in (b) 2-D and (c) 3-D simulations. The object is fixed at the centre of the domain. The ice has a circular cap at both ends in the 2-D simulations, and hemispherical caps for the 3-D cases.

3. Experimental and numerical methods

3.1. Experimental method and set-up

The experimental set-up consisted of an $l \times w \times d = {90}\,\textrm {cm} \times {50}\,\textrm {cm} \times {50}\,\textrm {cm}$ glass water tank, as sketched in figure 1(a). The tank was filled with water, and NaCl was dissolved to increase the salinity. Ambient water temperature and salinity were measured using a StarOddi DST CT probe with a temperature accuracy of ${\pm 0.1}\,\textrm {K}$ and a salinity accuracy of ${\pm 1}\,{\textrm {g}\,\textrm {kg}^{-1}}$ . The difference in composition between seawater and our NaCl solution is neglected in the salinity measurement, which adds an uncertainty of at most ${0.1}\,{\textrm {g}\,\textrm {kg}^{-1}}$ (Millero et al. Reference Millero, Feistel, Wright and McDougall2008). During all experiments, the ambient water temperature and salinity were measured to be within ${0.1}\,\textrm {K}$ and 0.1 $\textrm {g}\,\textrm {kg}^{-1}$ from the reported temperature and salinity, respectively. Prior to each experiment, we thoroughly mixed the water using a stirrer to achieve uniform temperature and salinity throughout the tank, then waited for the flow to quiesce. We used a 30 cm long ice cylinder, rounded on the bottom side, with a diameter of 5 cm for all experiments. Before each experiment, the ice was taken out of the freezer and left to equilibrate to room temperature for approximately 30 minutes, after which the ice core temperature was measured to be $T_s = {0.0}\,^\circ {\textrm {C}} \pm {0.5}\,\textrm {K}$ . The ambient conditions for all experiments are listed in table 1. Here, we used (2.8) to compute $\textit{Ra}_T$ and $\textit{Ra}_S$ , for which the density $\rho$ , kinematic viscosity $\nu$ and thermal diffusivity $\kappa _T$ were obtained from the empirical correlations of Sharqawy & Zubair (Reference Sharqawy and Zubair2010), evaluated at ambient temperature and salinity. Using the same expression for $\rho$ , the thermal expansion coefficient $\beta _T = - (1/\rho )\partial \rho /\partial T$ and the haline contraction coefficient $\beta _S = (1/\rho )\partial \rho /\partial S$ were obtained, and evaluated at ambient temperature and salinity.

Table 1.

Experimental conditions of the ambient water, including the ambient salinity $S_\infty$ and ambient temperature $T_\infty$ , and their corresponding Rayleigh numbers $\textit{Ra}_S$ and $\textit{Ra}_T$ , respectively. The last column also contains the density ratio $R_\rho$ . The dimensions of the cylinder were kept constant at $H = {30}\,\textrm {cm}$ and $D = {5}\,\textrm {cm}$ .

The shape evolution of the ice cylinder was quantified through boundary tracking, i.e. extraction of the ice–water interface from photographs taken at a regular time interval. This was done using a single camera pointing at the side of the cylinder, so only the 2-D projection of the cylinder was obtained. We used a Nikon D850 camera with a Sigma 105 mm macro objective, resulting in a resolution of ${42}\,{\unicode{x03BC} \textrm {m}\,\textrm {px}^{-1}}$ . We captured photographs at a 10 s interval for the experiments up to $\Delta S = {20}\,{\textrm {g}\,\textrm {kg}^{-1}}$ and a 5 s interval for the experiments with $\Delta S = {34}\,{\textrm {g}\,\textrm {kg}^{-1}}$ and $\Delta S = {81}\,{\textrm {g}\,\textrm {kg}^{-1}}$ , because of the higher melt rate. As illustrated in figure 1(a), the cylinder was lit from behind using a computer screen, on which the approximate shape of the cylinder was displayed. The width, height and curvature of this shape were manually adjusted during the experiment to ensure sufficient contrast along the edge of the cylinder for the entire duration of the experiment. Examples of photographs taken during the experiments are shown in figure 2. Each photograph was desaturated and then binarised using a global threshold, after which the edge contours were extracted using the border following algorithm developed by Suzuki & Abe (Reference Suzuki and Abe1985). Using these contours, the cylinder volume was computed as a vertical sum of disks, each with a height of 0.5 mm and a diameter based on the contours. The black regions in figure 2 show the presence of air bubbles inside the ice. We expect the effect of this entrapped air on our results to be minimal. By measuring the ice density, we estimated the volume fraction of air to be 1 %. Furthermore, most of the air is located in the centre of the cylinder, as the ice was frozen radially inward. In our analysis, we only use the 70 % of the volume that melts first, so we expect both the mean melt rate and the emergence of patterns to be unaffected by the entrapped air.

Figure 2.

Snapshots of nine experiments after 17 min of melting, with increasing ambient salinity from left to right. Black regions inside the ice are pockets of entrapped air.

3.2. Numerical method and set-up

Both 2-D and 3-D simulations are conducted in this study, with the ice represented as a cylinder in the 3-D simulations. For 2-D simulations, the flow is confined to a box of height $l$ and width $w$ , with the aspect ratio $\varGamma = w/l = 0.5$ . For 3-D simulations, the depth $d$ is the same as the width $w$ . No-slip, no-heat-flux and no-salt-flux boundary conditions are applied on all walls. A vertical ice cylinder with a height $H = 0.3l$ and diameter $D = 0.1w$ is fixed at the centre of the box and surrounded by uniformly warm and salty water, as shown in figure 1(b,c). The aspect ratios of the vertical ice cylinder and the box in simulations are consistent with those in experiments.

Simulations are performed using the second-order staggered finite difference code AFiD, which has been extensively validated and used to investigate a wide range of convection problems (Verzicco & Orlandi Reference Verzicco and Orlandi1996; Ostilla-Mónico et al. Reference Ostilla-Mónico, Yang, van der Poel, Lohse and Verzicco2015; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015; Yang, Verzicco & Lohse Reference Yang, Verzicco and Lohse2016). The extension of the AFiD code to study phase-change problems with the phase-field method has also been validated and discussed in recent studies (Liu et al. Reference Liu, Ng, Chong, Lohse and Verzicco2021; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023a , Reference Yang, Howland, Liu, Verzicco and Lohseb ).

In this study, several control parameters are fixed to restrict the extensive parameter space. Specifically, the Prandtl number and Schmidt number are set at $\textit{Pr}=1$ and ${Sc}=100$ , respectively, resulting in a Lewis number ${Le}=100$ . It is important to note that the Prandtl and Schmidt numbers ( $\textit{Pr}=1$ and ${Sc}=100$ ) used in the simulations are lower than those in the experiments, where $\textit{Pr}\approx 7$ and ${Sc}\approx 700$ . At higher Prandtl and Schmidt numbers, the thermal and saline boundary layers become significantly thinner, requiring substantially finer grid resolutions to accurately capture the boundary layer dynamics. Moreover, the mean melt rate decreases with increasing Prandtl and Schmidt numbers, leading to a longer melt time for melting a given fraction of the initial volume (or area) of the ice. Therefore, the choice of lower Prandtl and Schmidt numbers in the simulations is primarily motivated by the need to mitigate the substantial computational expense associated with resolving the thinner boundary layers and longer melt times characteristic of higher Prandtl and Schmidt numbers. In Appendix B, we examine the influence of varying Prandtl and Schmidt numbers in several cases with a low thermal Rayleigh number. While the mean melt rates are significantly lower at higher Prandtl and Schmidt numbers compared with lower values, the non-monotonic dependence of the mean melt rate on ambient salinity and the $\textit{Ra}_{S}$ scaling in the high-salinity regime remain unchanged across different Prandtl and Schmidt numbers. Thus, the use of lower Prandtl and Schmidt numbers in simulations does not alter the qualitative findings of the present study.

The initial temperature of the surrounding water is fixed at ${20}\,^\circ {\textrm {C}}$ for all simulations and the corresponding Stefan number is set to be ${Ste}=0.25$ . Our simulations cover a parameter range of ${1.5}\,\textrm {cm} \leqslant H\leqslant {30}\,\textrm {cm}$ and ${0}\,{\textrm {g}\,\textrm {kg}^{-1}}\leqslant \Delta S\leqslant {80}\,{\textrm {g}\,\textrm {kg}^{-1}}$ , roughly corresponding to $2.2\times 10^{5}\leqslant \textit{Ra}_{T}\leqslant 1.7\times 10^{9}$ and $0\leqslant \textit{Ra}_{S}\leqslant 2.3\times 10^{12}$ . The values of the height of cylinder $H$ and the corresponding $\textit{Ra}_T$ , as well as the range of $\Delta S$ and the corresponding $\textit{Ra}_S$ , are listed in table 2.

Table 2.

Control parameters used in the simulations, including the values of the height and diameter of cylinder $H$ and $D$ , and the corresponding $\textit{Ra}_T$ , the range of $\Delta S$ and the corresponding $\textit{Ra}_S$ as well as the density ratio $R_\rho$ . Note that $\textit{Ra}_S=0$ represents the melting in freshwater ( $\Delta S = {0}\,{\textrm {g}\,\textrm {kg}^{-1}}$ ). The notation $ [ a,\,b ]=\{x \in \mathbb{R}| a\leqslant x \leqslant b \}$ denotes a closed interval for $\Delta S$ and $\textit{Ra}_S$ .

Given that the diffusivity of salinity is significantly lower than that of temperature, the salinity field is resolved with much finer spatial resolution compared with that used for the temperature field. Additionally, the phase field is resolved with high precision to accurately capture the thin ice–water interface. Therefore, the multiple-resolution strategy of Ostilla-Mónico et al. (Reference Ostilla-Mónico, Yang, van der Poel, Lohse and Verzicco2015) is applied for both the salinity field and the phase field. For the 2-D cases with $H={1.5}\,\textrm {cm}$ and $\textit{Ra}_T=2.2\times 10^{5}$ , a uniform mesh of $n_{x}\times n_{y}=160\times 320$ is employed for the velocity and temperature fields. The uniform mesh for the salinity and phase fields is refined, increasing from $n_{x,f}\times n_{y,f}=480\times 960$ to $960\times 1920$ as $\Delta S$ rises from ${0}\,{\textrm {g}\,\textrm {kg}^{-1}}$ to ${80}\,{\textrm {g}\,\textrm {kg}^{-1}}$ . Here, $n_x$ and $n_y$ denote the horizontal and vertical grid resolutions for the velocity and temperature fields, respectively, while $n_{x,f}$ and $n_{y,f}$ represent the refined horizontal and vertical grid resolutions for the salinity and phase fields, respectively. When the height of the ice cylinder $H$ is increased to ${30}\,\textrm {cm}$ , a uniform mesh of $n_x \times n_y = 1920\times 3840$ is employed for the velocity and temperature fields. The refined uniform mesh for the salinity and phase fields ranges from $n_{x,f}\times n_{y,f}=5760\times 11\,520$ to $9600\times 19\,200$ as $\Delta S$ increases from ${0}\,{\textrm {g}\,\textrm {kg}^{-1}}$ to ${20}\,{\textrm {g}\,\textrm {kg}^{-1}}$ . Due to the high computational cost of the 3-D simulations, only seven 3-D cases were performed with $\textit{Ra}_{T}=2.2\times 10^{5}$ and $\Delta S={0, 1, 6, 20, 40, 60, 80}\,{\textrm {g}\,\textrm {kg}^{-1}}$ to confirm that the qualitative behaviours of the 3-D simulations are similar to those of the 2-D simulations. A uniform mesh of $n_x \times n_y \times n_z = 160 \times 320 \times 160$ is used for the velocity and temperature fields, and a finer uniform mesh ranging from $n_{x,f} \times n_{y,f} \times n_{z,f} = 480 \times 960 \times 480$ to $960 \times 1920 \times 960$ is employed for the salinity and phase fields, where $n_z$ and $n_{z,f}$ denote the grid resolution and refined grid resolution in the depth direction, respectively. The grid convergence test of the 2-D simulations is shown in Appendix A.

4. Melt rate

In this section, we examine the effect of the thermal and saline Rayleigh numbers on the mean melt rate. The normalised mean melt rates $\widetilde {f}$ as a function of $\textit{Ra}_{S_d}$ for different $\textit{Ra}_T$ in both experiments and simulations are shown in figure 3(a). The mean melt rate $\widetilde {f}$ exhibits a non-monotonic relation with $\textit{Ra}_{S_d}$ in both experiments and simulations: at constant $\textit{Ra}_{T}$ , as $\textit{Ra}_{S_d}$ increases, the mean melt rate $\widetilde {f}$ initially decreases and then increases again, which is similar to the observations by Yang et al. (Reference Yang, Howland, Liu, Verzicco and Lohse2023a ). This non-monotonic behaviour can be attributed to the competition between the thermal and saline buoyancy. At low ambient salinity, the saline buoyancy force is weak, resulting in a predominantly temperature-driven flow. As salinity increases, the saline buoyancy force strengthens and reaches equilibrium with thermal buoyancy at intermediate $\Delta S$ , leading to a reduced melt rate. Further increases in ambient salinity enhance the saline buoyancy force, leading to its dominance at high $\Delta S$ . Consequently, the meltwater moves upward more rapidly under the influence of strong saline buoyancy, accelerating the melt rate.

Figure 3.

(a) Normalised mean melt rate $\widetilde {f}$ as a function of $\textit{Ra}_{S_d}$ for different $\textit{Ra}_{T}$ in simulations as well as the normalised mean melt rate $\widetilde {f}$ as a function of $\textit{Ra}_{S}$ at $\textit{Ra}_{T}=5.3\times 10^{9}$ in experiments. (b) Relative mean melt rate $\widetilde {f}/\widetilde {f}_{0}$ as a function of $\textit{Ra}_{S_d}$ at $\textit{Ra}_{T}=2.2\times 10^{5}$ in 2-D and 3-D simulations, where $\widetilde {f}_{0}$ denotes the mean melt rate in freshwater ( $\Delta S={0}\,{\textrm {g}\,\textrm {kg}^{-1}}$ ) for the respective 2-D and 3-D cases. Note that $S_i$ is here assumed to be ${0}\,{\textrm {g}\,\textrm {kg}^{-1}}$ for the experiments.

At high $\textit{Ra}_{S_d}$ , while the saline boundary layer along the vertical cylinder remains laminar, the bulk flow becomes turbulent, where the mean salt flux follows $\overline {\textit{Nu}}_{S}\propto \textit{Ra}_{S_d}^{1/3}$ for the turbulent bulk flow with a laminar boundary layer (Bejan Reference Bejan1993; Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Holman Reference Holman2010). According to (2.11), we have $\widetilde {f} \propto \overline {\textit{Nu}}_{S}$ . Therefore, the mean melt rate $\widetilde {f}$ should follow $\widetilde {f}\propto \textit{Ra}_{S_d}^{1/3}$ at high $\textit{Ra}_{S_d}$ , which is in agreement with the results shown in figure 3(a). The deviation from this trend of the experimental data at high $\textit{Ra}_{S_d}$ is most likely due to the variation in temperature between experiments (table 1).

Additionally, some 3-D simulations of the vertical ice cylinder are conducted, as depicted in figure 3(a). It is found that the mean melt rates $\widetilde {f}$ in the 3-D simulations are larger than those in the 2-D simulations, primarily due to the larger value of the ratio between the surface area and the volume of the ice in the 3-D configuration. The mean melt rate $\widetilde {f}$ in the 3-D simulations exhibits a non-monotonic dependence on $\textit{Ra}_{S_d}$ , consistent with the trends observed in the 2-D simulations. At high $\textit{Ra}_{S_d}$ , the mean melt rate $\widetilde {f}$ appears to transition from a 1/4 scaling to a 1/3 scaling, rather than following a 1/3 scaling as observed in the 2D simulations. This suggests that the 1/3 scaling may emerge at a higher $\textit{Ra}_{S_d}$ in the 3-D configuration. Nevertheless, the similar non-monotonic behaviour observed in both 2-D and 3-D simulations indicates that the key physical mechanisms governing the melting process are well captured in the 2-D simulations.

The normalised mean melt rates $\widetilde {f}$ as a function of $\textit{Ra}_{T_d}$ for different $\Delta S$ as obtained from our simulations are shown in figure 4. We observe that the normalised mean melt rate $\widetilde {f}$ monotonically increases with increasing $\textit{Ra}_{T_d}$ for all $\Delta S$ considered. To gain further insight into the trend of the mean melt rate, the instantaneous vertical velocity $v_{y}$ , the temperature $T$ , salinity $S$ and fluctuating density $\rho ^{\prime }$ profiles in the horizontal direction at the mid-height of the vertical cylinder are shown in figure 5. The buoyancy forces induced by temperature and salinity act downward and upward, respectively, corresponding to a higher density at the wall compared with the ambient for temperature-driven flow (figure 5 b) and a lower density at the wall compared with the ambient for salinity-driven flow (figure 5 d).

Figure 4.

(a) Normalised mean melt rate $\widetilde {f}$ as a function of $\textit{Ra}_{T_d}$ for different $\Delta S$ in simulations. (b) Compensated normalised mean melt rate $\widetilde {f}{\textit{Ra}}_{T_d}^{-1/3}$ as a function of $\textit{Ra}_{T_d}$ for different $\Delta S$ .

When the ambient salinity is small, the temperature dependence dominates the buoyancy, causing the cold meltwater to descend, as indicated by the instantaneous vertical velocity $v_{y}$ in figure 5(a). This results in the formation of a laminar thermal boundary layer along the vertical ice cylinder within the laminar bulk flow, with the mean heat flux following the $\overline {\textit{Nu}}_{T}\propto \textit{Ra}_{T_d}^{1/4}$ scaling for the laminar bulk flow with a laminar boundary layer (Bejan Reference Bejan1993; Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Holman Reference Holman2010). According to (2.12), we have $\widetilde {f} \propto \overline {\textit{Nu}}_{T}$ . Consequently, the mean melt rate $\widetilde {f}$ adheres to the $\widetilde {f} \propto \textit{Ra}_{T_d}^{1/4}$ scaling relation for low $\Delta S$ values ( $\Delta S \lt {2}\,{\textrm {g}\,\textrm {kg}^{-1}}$ ), which is consistent with the observations seen in figure 4.

At intermediate $\Delta S$ , the buoyancy forces driven by temperature and salinity counteract, resulting in a bi-directional flow near the melt front, as depicted in figure 5(a). Due to the faster diffusion of heat compared with salt ( ${Le}=100$ ), the saline boundary layer is much thinner than the thermal boundary layer. Within the saline boundary layer, cold and fresh meltwater ascends due to salinity dominated buoyancy, while cold and saline ambient water descends due to thermal buoyancy outside the saline boundary layer and within the thermal boundary layer. This competition between thermal and saline buoyancies attenuates the buoyancy-driven flow and vertical velocity, consequently decelerating the melt rate. The flow around the vertical ice cylinder remains laminar, with the mean heat flux adhering to $\overline {\textit{Nu}}_{T}\propto \textit{Ra}_{T_d}^{1/4}$ (Bejan Reference Bejan1993; Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Holman Reference Holman2010). Accordingly, the mean melt rate $\widetilde {f}$ should also follow a $\widetilde {f}\propto \textit{Ra}_{T_d}^{1/4}$ scaling for intermediate $\Delta S$ ( ${2}\,{\textrm {g}\,\textrm {kg}^{-1}} \leqslant \Delta S \leqslant {10}\,{\textrm {g}\,\textrm {kg}^{-1}}$ ), consistent with the numerical findings presented in figure 4.

Figure 5.

(a) Normalised time-averaged vertical velocity $v_y/U_T$ as a function of the distance $x/D$ from the melt front at the mid-height of the vertical cylinder in simulations. The velocity is averaged over the time interval $0.2\leqslant t/t_{m}\leqslant 0.3$ . Here, $U_{T}=\sqrt {g\beta _{T}H\Delta T}$ is the free-fall velocity unit. Three cases with different ambient salinity ( $\Delta S={1}\,{\textrm {g}\,\textrm {kg}^{-1}}$ , ${3}\,{\textrm {g}\,\textrm {kg}^{-1}}$ , ${80}\,{\textrm {g}\,\textrm {kg}^{-1}}$ ) and $\textit{Ra}_T = 1.4\times 10^{7}$ are shown, with the same colours as in figure 4. (b–d) Temperature $T$ , salinity $S$ and fluctuating density $\rho ^{\prime }$ profiles as a function of the distance $x/D$ from the melt front at different $\Delta S$ .

However, when the ambient salinity is large, salinity dominates the buoyancy, causing the fresh meltwater to ascend, as depicted by the instantaneous vertical velocity $v_{y}$ in figure 5(a). Notably, figure 5(a) reveals that near the melt front ( $x/D\approx 0.06$ ), the magnitude of the positive $v_{y}$ in the case with $\Delta S={80}\,{\textrm {g}\,\textrm {kg}^{-1}}$ is nearly twice as large as those observed in the cases with $\Delta S={1}\,{\textrm {g}\,\textrm {kg}^{-1}}$ and ${3}\,{\textrm {g}\,\textrm {kg}^{-1}}$ . This increase can be attributed to the larger difference in density between the meltwater and the ambient water (figure 5 d), resulting in accelerated upward motion of the meltwater with increasing ambient salinity. The intrusion of warm ambient water onto the melt front occurs more rapidly, amplifying the melt rate and surface heat flux. With increasing $\Delta S$ , the surrounding boundary layer remains laminar, while the bulk flow begins to transition from laminar to turbulent, resulting in a transition in the scaling of the mean melt rate. Figure 4 shows this shift in scaling of the mean melt rate towards $\widetilde {f} \propto \textit{Ra}_{T_d}^{1/3}$ for high ambient salinity ( $\Delta S \gt {10}\,{\textrm {g}\,\textrm {kg}^{-1}}$ ), which is consistent with the mean heat flux scaling $\overline {\textit{Nu}}_{T}\propto \textit{Ra}_{T_d}^{1/3}$ for a turbulent bulk flow with a laminar boundary layer (Bejan Reference Bejan1993; Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Holman Reference Holman2010).

The $\textit{Ra}_S$ -scaling at high ambient salinity is analogous to the $\mathscr{V} \propto {Gr}^{1/3}=(\textit{Ra}/\textit{Pr})^{1/3}$ scaling derived by Mondal et al. (Reference Mondal, Gayen, Griffiths and Kerr2019), where $\mathscr{V}$ is the ablation velocity of a vertical ice face in cold seawater. In addition, the $\mathscr{V} \propto \Delta T_d^{4/3}$ relation for a vertical ice wall in cold seawater reported by Kerr & McConnochie (Reference Kerr and McConnochie2015) and Gayen et al. (Reference Gayen, Griffiths and Kerr2016) essentially represents a heat flux scaling of $\Delta T_d^{1/3}$ , similar to the $\textit{Ra}_{T_d}^{1/3}$ scaling that we observe here for high $\textit{Ra}_S$ . The similarity between the scalings reported by these studies for cold ambient water and those presented here for warm ambient water suggests that the physics in both cases is comparable, provided that the relative contributions of temperature and salinity on the buoyancy are the same.

5. Flow regimes

In this section, we investigate the influence of ambient salinity on the evolution of ice surface morphology. As discussed previously, the flow patterns can be divided into three regimes: temperature-driven, competing and salinity-driven flow. These flow regimes can be quantified using the density ratio $R_\rho$ , defined in (2.9). The limits $R_\rho \ll 1$ and $R_\rho \gg 1$ correspond to fully temperature-driven and fully salinity-driven flow, respectively. Using both experiments and simulations, we here explore the ice surface morphology dynamics for different $R_\rho$ by varying the ambient salinity while keeping the ambient temperature constant.

Figure 6.

Morphological evolution of the ice cylinder in simulations (left in each panel) and experiments (right in each panel) for the three different flow regimes: (a) temperature-driven flow; (b) competing flow; and (c) salinity-driven flow. The time intervals between the contours from simulations and experiments are $0.1 t_m$ and $0.15 t_m$ , respectively, where $t_m$ is the time needed to melt 70 % of the initial area in 2-D simulations and the volume in 3-D experiments. The innermost contours correspond to $t = t_m$ . The insets in panel (b) show an enlarged part of the contours, with the tracked position of the crests represented by red lines. In the experiments, the entire cylinder is not imaged and thus only the bottom part of each cylinder is visible.

Figure 7.

Snapshots of the temperature (left) and salinity (right) fields for the three different flow regimes in the simulations at time $t/t_{m}=0.8$ : (a) temperature-driven flow, (b) competing flow, and (c) salinity-driven flow. Here $\widetilde {T}=T/\Delta T$ and $\widetilde {S}=S/\Delta S$ .

Both our experiments and simulations show that the evolution of the ice surface morphology depends heavily on $R_\rho$ , exhibiting distinct dynamics for each of the three flow regimes (figure 6). Remarkably, the dynamics of the morphology from the 2-D simulations closely resembles those obtained from 3-D experiments at similar $R_\rho$ . In the following paragraphs, we discuss the evolution of the ice morphology for each flow regime separately and relate it to the corresponding flow, based on the instantaneous temperature and salinity fields obtained from simulations (figure 7).

For temperature-driven flow ( $R_\rho \lt 1$ ), we found that the ice cylinder exhibits a pointed top and a flat bottom (figure 6 a). The denser cold meltwater adjacent to the cylinder descends along its surface, insulating the lower part of the cylinder while exposing the upper part to the warmer ambient water. This results in a higher melt rate at the top as compared with the bottom, leading to the observed morphological dynamics. The difference in melt rate between the top and the bottom is obvious in the contours from the simulation (left in figure 6 a), but less so for the experiment (right in figure 6 a). This is possibly due to the difference in shape of the top of the cylinder, which is round in the simulation, but flat in the experiment (for practical reasons). As this will only affect a downward flow, it will mainly be relevant in the temperature-driven regime. Additionally, a cavity is formed at the bottom of the ice cylinder in the simulation (left in figure 6 a for times $0.1 \leqslant t/t_m \leqslant 0.8$ ). This cavity formation can be attributed to flow separation at the bottom of the cylinder, which generates attached convective rolls. These convection rolls, evident from the instantaneous salinity field (figure 7 a), enhance local mixing, increasing the local heat flux and therewith the melt rate, ultimately leading to the formation of the cavity. Around $t = 0.8 t_m$ , the diameter of the cylinder is almost half its initial value, such that the convection rolls are disturbed by the downward flow along the cylinder, causing them, and therefore the cavity, to disappear. In the experimental contours (right in figure 6 a), the cavity is not visible due to the optical limitations of the side view.

For competing flow ( $R_\rho \approx 2$ ), a regular scallop pattern appears on the side of the cylinder (figure 6 b). Between the downward-flowing layer of cold ambient water and the ice surface, a thin layer of fresh meltwater flows upward along the side of the cylinder. This bi-directional flow is clearly evident from the instantaneous vertical velocity profile for $\Delta S={3}\,{\textrm {g}\,\textrm {kg}^{-1}}$ (figure 5 a). The shear between the upward- and downward-flowing layers induces a Kelvin–Helmholtz-type instability, akin to that observed by Weady et al. (Reference Weady, Tong, Zidovska and Ristroph2022) when they melted an ice cylinder in ${5.6}\,^\circ {\textrm {C}}$ freshwater. Similar to the case of Weady et al. (Reference Weady, Tong, Zidovska and Ristroph2022), our results show that the scallop pattern is carved by an unstable bidirectional flow. In their case, the flow consists of a positively buoyant inner layer of water below the maximum density temperature, nested inside a negatively buoyant outer layer of water above the maximum density temperature. In our case, the positively buoyant inner layer is cold and fresh meltwater, and the negatively buoyant outer layer is cold and saline ambient water. Due the difference in salinity between the layers, the density difference is 10–100 times larger than in the freshwater case of Weady et al. (Reference Weady, Tong, Zidovska and Ristroph2022). In addition, the widths of the inner and outer layers are controlled by the saline and thermal diffusivity, respectively, while only the latter is relevant in freshwater. Despite these differences, the resulting scallop pattern observed in our saline water cases is qualitatively similar to the freshwater case of Weady et al. (Reference Weady, Tong, Zidovska and Ristroph2022). We conclude that scallop patterns emerge due to the presence of a bidirectional flow, regardless of the nature of this flow. The characteristics of these scallops are discussed in more detail in § 6.

In salinity-driven flow ( $R_\rho \gg 1$ ), the ice melts fastest just above the bottom of the cylinder, narrowing the cylinder’s width at this point (figure 6 c). Above this point, cold meltwater ascends along the cylinder, shielding the ice from warm ambient water and consequently reducing its melting rate. Below this point, thermal driving is strong enough to prevent meltwater from flowing up along the cylinder. The cold and fresh meltwater ascends within the saline boundary layer, while the cold and saline ambient water descends outside the saline boundary layer but within the thermal boundary layer, as confirmed by the instantaneous snapshots of temperature and salinity fields (figure 7 c). This bi-directional flow shields the melt front from direct contact with the warm ambient water, leading to a slower local melt rate at the bottom region. At the point of minimum width, the melt rate increases due to the intrusion of warm ambient water, as cold meltwater either flows up along the cylinder or is advected downward by the temperature-driven flow. This observation is similar to the bifurcation point observed by Josberger & Martin (Reference Josberger and Martin1981) in their experiments on a vertical ice wall. Furthermore, similar to the cavity formed in the cylinder with temperature-driven flow (figure 6 a), a cavity is formed with salinity-driven flow (figure 6 c), but at the top of the cylinder, instead of at the bottom. Here, the upward flow of meltwater separates from the ice surface at the top, generating convection rolls that enhance heat transport to the ice and thus the melt rate.

The extent of these three flow regimes and their implications on the ice morphology are summarised in the $\textit{Ra}_{T}$ – $\textit{Ra}_{S}$ and $\textit{Ra}_{T}$ – $R_{\rho }$ phase diagrams in figures 8(a) and 8(b), respectively. The boundary between the temperature-driven and competing regimes is defined at $R_\rho = 1$ , corresponding to the zero wall-shear condition established by Carey & Gebhart (Reference Carey and Gebhart1982). For the boundary between the competing and salinity-driven regimes, we adopt $R_\rho = 6.37$ , corresponding to the density ratio for which Josberger & Martin (Reference Josberger and Martin1981) found the wall-normal integral of buoyancy near a vertical ice wall to be zero. It should be noted that these values of $R_\rho$ serve as indicative rather than strict boundaries, as the transitions between the regimes are continuous. Nevertheless, our observations of scallops in both experiments and simulations align closely with the competing flow regime (figure 8). Scallop patterns are only observed for $\textit{Ra}_T \gtrsim 10^7$ , indicating a minimum thermal driving for triggering an instability. Bottom and top cavities only appear at sufficiently high $\textit{Ra}_T$ and $\textit{Ra}_S$ , respectively, which is likely linked to a critical Reynolds number and thus a minimal Rayleigh number for which flow separation occurs.

Figure 8.

(a) $\textit{Ra}_{T}$ – $\textit{Ra}_{S}$ and (b) $\textit{Ra}_{T}$ – $R_{\rho }$ phase diagrams for various morphological patterns observed in simulations and experiments. The green, orange and purple regions correspond to the temperature-driven, competing and salinity-driven flow regimes, respectively. The large hexagonal markers in panel (b) indicate the minimum mean melt rate case for each $\textit{Ra}_T$ , as discussed in § 4.

6. Properties of the ice scallops

In this section, we discuss the scallop patterns observed in the thermal-saline competing regime (figures 6 b and 7 b). The typical length scale for the scallop patterns, wavelength $\lambda$ , obviously depends on the ambient conditions. To investigate the dependence of the scallop wavelength on the ambient conditions, we compute the mean (in space and time) of all scallop wavelengths $\bar {\lambda }$ for each case. The mean wavelengths $\bar {\lambda }$ as a function of the density ratio $R_{\rho }$ in simulations and experiments are depicted in figure 9(a). As shown by the data from simulations, the mean wavelength $\bar {\lambda }$ increases with increasing $\textit{Ra}_T$ at similar $R_\rho$ , indicating that the thermal driving strength influences the scallop pattern. At similar $\textit{Ra}_T$ , the wavelength observed at low $R_\rho$ (low salinity) is up to twice the value observed for high $R_\rho$ (high salinity). Yet, no clear correlation between $\bar {\lambda }$ and $R_\rho$ emerges, given that the standard deviation spread in wavelengths exceeds 10 % for each case, as indicated by the interval markers in figure 9(a). The mean wavelengths obtained from the simulations with $\textit{Ra}_T = 1.7 \times 10^9$ agree well with those observed in experiments performed at similar ambient conditions ( $\textit{Ra}_T = 5.3 \times 10^9$ ).

Figure 9.

(a) Mean wavelengths $\bar {\lambda }$ , (b) compensated mean wavelengths $ ( \bar {\lambda } /H )\textit{Ra}_{T}^{1/6}$ , and (c) scallop migration velocities $\chi$ as a function of the density ratio $R_{\rho }$ in simulations and experiments. Each point corresponds to a single simulation or experiment. The interval markers represent one standard deviation spread in wavelength or migration velocity. The inset in panel (a) shows the definition of the wavelength $\lambda$ and the coordinate system for computation of $\chi$ in (6.3).

To further explore the dependence of $\bar {\lambda }$ on $\textit{Ra}_T$ and $R_\rho$ , we perform a least-squares fit of the form $\bar {\lambda }/H = a \textit{Ra}_T^b R_\rho ^c$ , where $H$ is the cylinder height, and $a$ , $b$ and $c$ are the free parameters. It yields the empirical relation:

(6.1) \begin{equation} \frac {\bar {\lambda }}{H} = (3.8 \pm 1.8) \; \textit{Ra}_T^{-0.16 \,\pm\, 0.02}\, \; R_\rho ^{-0.54 \,\pm\, 0.15}, \end{equation}

where the $\pm$ -values correspond to standard errors.

A possible explanation for this scaling of $\bar {\lambda }$ with $\textit{Ra}_T$ can be given by constructing a typical length scale for the temperature-driven flow. For a vertical cylinder in a single-phase fluid with $\textit{Pr} \geqslant 1$ , the velocity $U$ of a thermally driven flow follows $U \propto (\kappa _T/H) \textit{Ra}_T^{1/2}$ (Bejan Reference Bejan1993), with the local buoyancy time scale $\tau = (\nu /(g\beta _T\Delta T)^2)^{1/3}$ . Assuming that the mean wavelength is proportional to the typical length scale of the flow $U\tau$ , we obtain

(6.2) \begin{equation} \bar {\lambda } \propto U \tau = \frac {\kappa _T}{H} \textit{Ra}_T^{1/2} \left (\frac {\nu }{(g\beta _T\Delta T)^2}\right )^{1/3} = H\textit{Ra}_T^{-1/6}. \end{equation}

The scaling in (6.2) can be evaluated using the compensated mean wavelengths $ ( \bar {\lambda } /H )\textit{Ra}_{T}^{1/6}$ , as shown in figure 9(b). We find that this provides a better collapse of the mean wavelengths for the parameter space explored in this work. Moreover, the exponent obtained from the fit ( $\textit{Ra}_T^{-0.16 \,\pm\, 0.02}$ ) is in good agreement with the scaling $\bar {\lambda }/H \propto \textit{Ra}_T^{-1/6}$ in (6.2). Given the large spread in wavelengths, this scaling clearly warrants further investigation. Specifically, the differences and interactions between individual scallops should be explored. Scallop wavelengths could depend on their vertical position (Weady et al. Reference Weady, Tong, Zidovska and Ristroph2022), and neighbouring scallops could merge or split over time. Such details are not accessible when a mean wavelength is used. Although a robust theory for the pattern formation therefore remains elusive, we infer that both the driving strength and the ratio of thermal to saline driving play important roles in determining the scallop wavelength, underscoring their significance in future investigations.

Another notable feature of the scallop pattern is that it migrates downward over time (figure 6 b). Similar migration of scallops has been observed in horizontal blocks of dissolving plaster and melting ice, and inclined blocks of dissolving caramel, salt and plaster, subject to forced convection (Blumberg & Curl Reference Blumberg and Curl1974; Gilpin et al. Reference Gilpin, Hirata and Cheng1980; Bushuk et al. Reference Bushuk, Holland, Stanton, Stern and Gray2019; Cohen et al. Reference Cohen, Berhanu, Derr and Courrech du Pont2020).

To quantify the downward migration, we define a dimensionless migration velocity $\chi$ as the vertical movement of the scallop pattern compared with its horizontal movement. The vertical speed is given by the change in vertical position $y_c(t)$ of scallop crests (maxima in local radius) over time ( ${\textrm d} y_c/{\textrm d} t$ ), and the horizontal speed is given by the height-averaged horizontal melt rate ( $\left \langle \partial r/\partial t \right \rangle _{y}$ ), leading to the following definition of the migration velocity $\chi$ :

(6.3) \begin{equation} \chi =\frac {\left \langle \frac {{\textrm d}{y}_c}{{\textrm d} t}\right \rangle _{t, N}}{\left \langle \frac {\partial r(y, t)}{\partial t} \right \rangle _{y, t}}, \end{equation}

Figure 10.

(a) Normalised horizontal melt rate along the contour of the cylinder from the experiment at $R_\rho = 1.4$ ( $\textit{Ra}_T = 4.2\times 10^{9}$ , $\textit{Ra}_S = 6.0 \times 10^{11}$ ), averaged over 10 min. Positive and negative values indicate above average and below average melt rates, respectively. (b) Snapshot of the temperature field from the numerical simulation at $R_\rho = 2.0$ ( $\textit{Ra}_T = 1.7 \times 10^9$ , $\textit{Ra}_S = 3.4 \times 10^{11}$ ), with local melt rates on the contour, averaged over 2 minutes. Points A and B indicate the accumulation of cold meltwater and intrusion of warm ambient water, respectively, also sketched in the inset of panel (a).

where the subscripts $y$ , $t$ and $N$ indicate averages over height, time and number of scallop crests, respectively. The coordinate system and crest locations are sketched in the inset of figure 9(a). Blumberg & Curl (Reference Blumberg and Curl1974) and Bushuk et al. (Reference Bushuk, Holland, Stanton, Stern and Gray2019) observed a similar pattern migration in their experiments and quantified this movement using a crest propagation angle $\phi$ . The migration velocity $\chi$ (6.3) is related to this propagation angle as $\chi = \cot (\phi )$ . With increasing $R_\rho$ , the migration velocity decreases monotonically (figure 9 c). The average migration velocities from simulations agree well with those observed in experiments, further supporting the suitability of the 2-D simulations in representing the 3-D experiments. There is little variation between different $\textit{Ra}_T$ , and thus, $\chi$ is nearly independent of the thermal driving strength.

To explain the observed migration, we show the normalised local melt rate along the ice contour for the experiment at $R_\rho = 1.4$ in figure 10(a). It is defined as $\hat {m} = (\langle \dot {r} \rangle _t - \langle \dot {r} \rangle _{y, t} ) / \langle \dot {r} \rangle _{y, t}$ , where the subscripts $y$ and $t$ indicate averages over height and time, respectively. For the computation of $\langle \dot {r} \rangle _t$ , the horizontal melt rate profile at each time $t$ was shifted vertically by a distance $-\langle {\textrm d}y_c/{\textrm d}t\rangle _{t,N} (t-t_0)$ to account for the scallop migration since a time $t_0$ . The melt rates $\hat {m}$ shown in figure 10 are therefore relative to the scallop position. As shown in figure 10(a), the melt rate is larger at the upper side of each crest compared with its lower side. This asymmetry leads to faster ablation of the ice on the upper side of the crests, resulting in the apparent downward migration of the pattern.

To explore the origin of this asymmetric melt rate, we show the instantaneous temperature field from the simulation with $R_\rho = 2.0$ in figure 10(b). A notable temperature disparity is evident between the two sides of the crests, resulting in an uneven melt rate. Between adjacent crests, cold meltwater moves up along the ice and accumulates on the lower side of a crest (points A in figure 10 a,b). This displacement of fluid allows warm ambient water to intrude (points B in figure 10 a,b).

Using this proposed mechanism for pattern migration, we can explain our observations from figure 9(c). With increasing $R_\rho$ , the saline buoyancy of the meltwater increases, such that less meltwater accumulates behind the scallop crests and instead moves upward past the crests. This reduces the local differences in melt rate and therefore decreases the migration velocity, explaining the decrease of $\chi$ with increasing $R_\rho$ in figure 9(c).