## 1. Introduction

Melting and freezing have huge relevance in various fields, with a wide range of applications in nature and technology, including sea ice (Holland, Bitz & Tremblay Reference Holland, Bitz and Tremblay2006), phase-change materials (Dhaidan & Khodadadi Reference Dhaidan and Khodadadi2015), aircraft icing (Cao, Tan & Wu Reference Cao, Tan and Wu2018), icebergs (Ristroph Reference Ristroph2018) and icy moons (Spencer *et al.* Reference Spencer, Pearl, Segura, Flasar, Mamoutkine, Romani, Buratti, Hendrix, Spilker and Lopes2006; Kang & Flierl Reference Kang and Flierl2020). Accurately quantifying the melt rate of glacial ice in the ocean is vital for constraining estimates of sea level rise under various climate change scenarios (Edwards *et al.* Reference Edwards2021; Cenedese & Straneo Reference Cenedese and Straneo2023). The presence of salt in terrestrial (and possibly extraterrestrial) oceans introduces double-diffusive convection driven by both temperature and salinity variations (Kang *et al.* Reference Kang, Mittal, Bire, Campin and Marshall2022). The coupling of such flows to a moving phase boundary from a mathematical point of view constitutes a so-called Stefan problem (Rubinstein Reference Rubinstein1971).

To better understand such highly complex systems, we consider a sufficiently simplified model problem that still contains the rich phenomenology of the turbulent flow around the melting ice observed in reality. Much previous work has focused on melting with *single*-component convective flows, where the melting dynamics and convection are determined solely by temperature variations (Davis, Müller & Dietsche Reference Davis, Müller and Dietsche1984; Favier, Purseed & Duchemin Reference Favier, Purseed and Duchemin2019; Yang *et al.* Reference Yang, Howland, Liu, Verzicco and Lohse2023). Extensions to this approach, using both experiments and simulations, have also considered the effects of shear on morphology and temperature/salt (Couston *et al.* Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2021; Hester *et al.* Reference Hester, McConnochie, Cenedese, Couston and Vasil2021) and the effects of rotation (Ravichandran & Wettlaufer Reference Ravichandran and Wettlaufer2021; Toppaladoddi Reference Toppaladoddi2021; Ravichandran, Toppaladoddi & Wettlaufer Reference Ravichandran, Toppaladoddi and Wettlaufer2022) on the phase transition process, as well as the dependence on the initial conditions (Purseed *et al.* Reference Purseed, Favier, Duchemin and Hester2020). The pioneering work of Wang, Calzavarini & Sun (Reference Wang, Calzavarini and Sun2021*a*) and Wang *et al.* (Reference Wang, Calzavarini, Sun and Toschi2021*b*,Reference Wang, Jiang, Du, Sun and Calzavarini*c*) on the density maximum effect has identified distinct regimes of the flow pattern, ice morphology and melt speed under different temperatures and geometries.

However, salinity significantly complicates the problem, as it modifies the density and the melting point of aqueous ice. The importance of salinity in ice melting has been demonstrated experimentally by Huppert & Turner (Reference Huppert and Turner1978, Reference Huppert and Turner1980) and McConnochie & Kerr (Reference McConnochie and Kerr2016). They showed that the meltwater spreads into the liquid in a series of horizontal layers. Also, the ice forms layered structures corresponding to the flow structures, and even more complex mushy layer structures can form (Worster Reference Worster1997), where convection also occurs (Du *et al.* Reference Du, Wang, Jiang, Calzavarini and Sun2023). Numerical simulations have been used to study the layer structures in laterally cooled double-diffusive convection, with a temperature gradient in the horizontal direction and a salinity gradient in the vertical direction (Kranenborg & Dijkstra Reference Kranenborg and Dijkstra1998; Chong *et al.* Reference Chong, Yang, Wang, Verzicco and Lohse2020). However, the coupling of such a flow and the melting process have not been numerically modelled, owing to the challenge of properly representing the effect of salinity on ice melting, as well as computational limitations. Here we will overcome these limitations to quantitatively answer how salinity affects the melt rate and shape evolution of ice.

We conduct numerical simulations of a fixed vertical ice block, melting from the side by salt-stratified water. The Navier–Stokes equations and the advection equations for temperature and salinity are coupled to the phase-field for the ice–water interface, a model which is widely used for phase boundary evolution (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, Chong, Liu, Verzicco and Lohse2022*a*). Layered structures on the melt front are observed and quantitatively agree with the experiments of Huppert & Turner (Reference Huppert and Turner1980). Furthermore, a non-monotonic trend is observed, where the melt rate is first reduced and then enhanced as the salinity in the water increases. Despite the complexity of the interaction between the moving boundary and the turbulent flow, we provide a simple theoretical model based on a balance of buoyancy forces, which predicts the dependence of the minimal melt rate on salinity and temperature throughout.

## 2. Numerical method and set-up

We implement the phase-field method presented by Hester *et al.* (Reference Hester, Couston, Favier, Burns and Vasil2020) to simulate the melting solid in a multi-component fluid, in which 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 location of the interface defined implicitly by $\phi =1/2$. The governing equations, constrained by the incompressibility condition $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} = 0$, read

Here, $\boldsymbol {u}$ is the velocity vector, $p$ is the kinematic pressure, $g$ is gravitational acceleration in the vertical direction $\boldsymbol {e_z}$, $\nu$ is the kinematic viscosity of the liquid, $\kappa _T$ and $\kappa _S$ are the diffusivities of heat and salt respectively, $\mathcal {L}$ is the latent heat, and $c_p$ is the specific heat capacity. The limit $\varepsilon \to 0$ gives the exact boundary conditions for $T$ and $S$:

The liquidus slope $m=0.056\,^\circ {\rm C} ({\rm g}\ {\rm kg}^{-1})^{-1}$ describes how the melting point depends on salinity. Here, we neglect the density difference between the two phases, which would lead to a Stefan flow, that is, a non-zero fluid velocity at the interface. Based on the melt rate and flow velocity magnitude from our simulations, we estimate this flow to be multiple orders of magnitude weaker than the buoyancy-driven flow near the interface.

The flow is confined to a rectangular box of height $H$ and aspect ratio $\varGamma =L_x/H=1$. For the three-dimensional (3D) simulations, the depth-wise aspect ratio is set to $\varGamma _y = L_y/H = 0.5$. No-heat-flux, no-salt-flux and no-slip boundary conditions are applied on all walls. Initially, we place a vertical ice block with a thickness of $0.1H$ and uniform temperature at the melting point, as shown in figure 2(*a*). The initial temperature field $T$ and salinity field $S$ are prescribed as follows, with a linear salinity profile in the vertical ($z$) direction for the liquid phase:

*a*,

*b*)\begin{equation} T= \begin{cases} T_w, & x<0.9H, \\ T_i, & x\ge 0.9H, \end{cases} \quad S = \begin{cases} S_{bot} + (S_{top} - S_{bot}) z/H, & x<0.9H, \\ 0, & x\ge 0.9H. \end{cases} \end{equation}

The initial solid temperature $T_i = 0\,^\circ$C is set to the equilibrium melting temperature for zero salinity. $S_{top}$ and $S_{bot} \ge S_{top}$ are the initial values of salinity at the top and bottom boundaries, respectively. Since the initial linear background profile is incompatible with the no-salt-flux boundary, thin diffusive boundary layers with zero gradient in salinity form at the boundaries. However, the slow diffusion of salt means that these boundary layers remain thin throughout and hardly affect the dynamics. From these initial conditions, we can define a temperature scale and two salinity scales, accounting for variations in the vertical due to stratification, and in the horizontal between the ice (which has $S_i=0$) and the mean salt concentration in the water:

*a*–

*c*)\begin{equation} {\rm \Delta} T = T_w - T_i, \quad {\rm \Delta} S_{v} = S_{top} - S_{bot}, \quad S_m = (S_{top} + S_{bot})/2 . \end{equation}

Based on the commonly used Oberbeck–Boussinesq approximation, which retains the density variation only in the buoyancy term, we use a simplified, yet realistic, equation of state (EOS) for water at atmospheric pressure (Roquet *et al.* Reference Roquet, Madec, Brodeau and Nycander2015), defined as

where $\rho '=\rho - \rho _0$ is the fluid density perturbation from a reference value $\rho _0$, and the coefficients have values $C_b=0.011\ {\rm kg}\ {\rm m}^{-3}\ {\rm K}^{-2}$, $b_0=0.77\ {\rm kg}\ {\rm m}^{-3} ({\rm g}\ {\rm kg}^{-1})^{-1}$, $T_0=4\,^\circ$C, $c_S=-0.25\ {\rm K}\ ({\rm g}\ {\rm kg}^{-1})^{-1}$. The temperature-dependence of density perturbations arising from (2.8) is shown in figure 2(*b*) for a range of fixed salinity values. We shall study the initial value problem of how the ice melts (interface pattern and overall melt rate) in different salinity environments.

The applied phase-field model was initially derived based on entropy conservation, which guarantees thermodynamic consistency, and also satisfies the Gibbs–Thomson relation (Wang *et al.* Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993; Favier *et al.* Reference Favier, Purseed and Duchemin2019; Hester *et al.* Reference Hester, Couston, Favier, Burns and Vasil2020). Following Hester *et al.* (Reference Hester, Couston, Favier, Burns and Vasil2020), the related numerical parameters include $\varepsilon$, which is used to measure the diffusive interface thickness and set to be equal to the grid spacing, and the penalty parameter $\eta$, which is used to damp the velocity to zero in the solid phase and is set equal to the time step. To ensure the accuracy of the simulation, we further impose a direct forcing of the velocity to zero for $\phi >0.9$, preventing spurious motions inside the solid phase (Howland Reference Howland2022). The diffusivity of the phase field is defined by $D=6\gamma \kappa _T/5\varepsilon \mathcal {L}$, following Hester *et al.* (Reference Hester, Couston, Favier, Burns and Vasil2020), where $\gamma$ is a surface energy coefficient related to the Gibbs–Thomson effect. We set $D=1.2\kappa$ and $\varepsilon {\rm \Delta} T/\gamma = 1.0$. A more detailed discussion of the parameters and the validation of the phase-field method can be found in our own software documentation (Howland Reference Howland2022) and in Hester *et al.* (Reference Hester, Couston, Favier, Burns and Vasil2020).

With the parameters of the phase-field model fixed, the system is uniquely determined by six dimensionless control parameters. For these we take the thermal and solutal Rayleigh numbers, the Prandtl number, the Schmidt number, the Stefan number and the density ratio between vertical and horizontal salinity difference $\varLambda _S$:

*a*–

*c*)\begin{gather} Ra_T = \frac{g C_b {\rm \Delta} T^2 H^3}{2\nu\kappa_T}, \quad Ra_S = \frac{g b_0 S_m H^3}{ \nu\kappa_S}, \quad Pr = \frac{\nu}{\kappa_T}, \end{gather}

*a*–

*c*)\begin{gather}Sc = \frac{\nu}{\kappa_S}, \quad St = \frac{\mathcal{L}}{c_p {\rm \Delta} T}, \quad \varLambda_S =\frac{{\rm \Delta} S_v}{S_m} . \end{gather}

Furthermore, we can define the Lewis number as the ratio of heat and salt diffusivity, as well as the density ratio between temperature and horizontal salinity difference $\varLambda _T$:

*a*,

*b*)\begin{equation} Le = \frac{\kappa_T}{\kappa_S} = \frac{Sc}{Pr}, \quad \varLambda_T =\frac{2b_0S_m}{C_b{\rm \Delta} T^2} = \frac{Ra_S }{Ra_T}. \end{equation}

Because the parameter space is so large, some of the control parameters have to be fixed in order to make the study feasible. We fix $Pr=10$ and $Sc=1000$ (i.e. $Le=100$) as relevant values for seawater in all cases. Our simulations cover a parameter range of $10\,^\circ {\rm C}\le {\rm \Delta} T\le 20\,^\circ {\rm C}$, $0\le S_m\le 15\ {\rm g} \ {\rm kg}^{-1}$ and $2.5\ {\rm cm}\le H\le 10\ {\rm cm}$, corresponding roughly to $10^6 \le Ra_T\le 10^8$, $0 \le Ra_S\le 2.5\times 10^{10}$ and $4\le St\le 8$. Unless otherwise specified, we fix the initial temperature of the water as $T=20\,^\circ$C and the domain height as $H=5\ {\rm cm}$ (corresponding to $Ra_T=10^7$).

Simulations are performed using the second-order staggered finite difference code AFiD (van der Poel *et al.* Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015), which has been extensively validated (Kooij *et al.* Reference Kooij, Botchev, Frederix, Geurts, Horn, Lohse, van der Poel, Shishkina, Stevens and Verzicco2018) and used to study a wide range of convection problems, including double-diffusive convection (Yang, Verzicco & Lohse Reference Yang, Verzicco and Lohse2016; Yang *et al.* Reference Yang, Verzicco, Lohse and Caulfield2022*b*). The extension of the AFiD code to include two phases approached with the phase-field method was discussed and validated in Liu *et al.* (Reference Liu, Ng, Chong, Lohse and Verzicco2021), Yang *et al.* (Reference Yang, Chong, Liu, Verzicco and Lohse2022*a*).

Since salt diffuses much more slowly than heat, the resolution for the salinity is more demanding than that for the temperature, which limits the simulations to moderate $Ra_T$ and $Ra_S$. Here, the multiple-resolution strategy of Ostilla-Mónico *et al.* (Reference Ostilla-Mónico, Yang, Van Der Poel, Lohse and Verzicco2015) is employed for both salinity and phase field, since it is vital for the accuracy of the phase-field method that the interface thickness be smaller than the thin salt concentration boundary layer. Specifically, for the case of $Ra_T=10^8$ and $Ra_S=2\times 10^{10}$, we use a uniform mesh of $432\times 432$ for the velocity and temperature field and $2880\times 2880$ for the salinity field and phase field. For the 3D case of $Ra_T=1.5\times 10^7$ and $Ra_S=2\times 10^{9}$, we use a uniform mesh of $288\times 288\times 144$ and a refined mesh of $864\times 864\times 432$. A grid-independence test has been conducted to ensure that the same conclusions are obtained, as shown in figure 1.

## 3. Effect of salinity on the structure of the melting interface

To reveal the effect of salinity on the shape evolution and the melt rate of the ice, we begin with a qualitative description of how the melt front shape depends on the vertical salinity gradient ${\rm \Delta} S_v$. In figure 2(*c*–*e*) typical temperature and salinity fields for three different 3D simulations are shown, where we fix ${\rm \Delta} T=20\ {\rm K}$, $H=5\ {\rm cm}$, $S_m=5\ {\rm g} \ {\rm kg}^{-1}$, and vary ${\rm \Delta} S_v$.

At relatively low or zero vertical salinity variation ${\rm \Delta} S_v =0$ (figure 2*c*), when there is no initial stable stratification, salinity and temperature are mostly uniform in the bulk, and a concave ice melt front forms because of the large-scale circulation, as shown in figure 2(*c*), right. At moderate ${\rm \Delta} S_v (=5\ {\rm g}\ {\rm kg}^{-1})$ (figure 2*d*), a layered structure occurs for the temperature and salinity fields, and correspondingly also for the ice melt front, similarly to the experiments of Huppert & Turner (Reference Huppert and Turner1980). With a further increase in ${\rm \Delta} S_v(=10\ {\rm g}\ {\rm kg}^{-1})$ (figure 2*e*), more layers appear in the liquid phase as compared to the case of ${\rm \Delta} S_v=5\ {\rm g}\ {\rm kg}^{-1}$, while the layer structure disappears on the melt front, since the convective flow at the interface is weakened by the stronger stable stratification. In figure 2(*e*), the top salinity value approaches zero, causing the fresh meltwater plume to intrude at the base of the broad upper layer. This colder plume, with the same salinity as the layer, enters at a level of neutral buoyancy. By contrast, for figure 2(*c*,*d*), the less saline meltwater forms a thin layer at the top of the domain. The density contrast between this thin upper layer and the ambient creates strong local stratification, hindering circulation and heat transport to the upper ice. This results in a concave shape at the top. In figure 2(*e*), warm water is transported to the upper ice through convection.

With this qualitative behaviour established, we now focus on a more detailed analysis of how both vertical and horizontal variations in salinity affect the system. In figure 3(*a*), we present the temperature and salinity fields from two-dimensional (2D) simulations at various horizontal and vertical salinity variations $S_m$ and ${\rm \Delta} S_v$. One can see that the structure of the flow and melt front depends mainly on ${\rm \Delta} S_v$. Therefore, the vertical salinity gradient and horizontal temperature gradients are the main driving factors of the flow structure, consistent with the findings of Huppert & Turner (Reference Huppert and Turner1980) and Chong *et al.* (Reference Chong, Yang, Wang, Verzicco and Lohse2020). Note that Weady *et al.* (Reference Weady, Tong, Zidovska and Ristroph2022) found a similar layer pattern in ice melting in fresh water near $4\,^\circ$C, which is explained by the density anomaly for fresh water. Unlike in fresh water, the layers in salty water are due to the competition between thermal and salinity-driven buoyancy.

The different diffusivities of heat and salt play a significant role in the dynamics outside the thin fresh layer. Since heat diffuses much faster than salt ($Le=100$), a region of cold, saline water is produced, which sinks because of thermal-induced buoyancy. At ${\rm \Delta} S_v = 0$, the lack of salt stratification allows this water to sink. The accumulation of cold water causes thicker ice at the bottom than in the middle; see the leftmost column in figure 3(*a*). Combined with the accumulation of fresh water at the top, which weakens the flow, this results in a concave shape for the ice, with the thinnest ice in the middle.

At ${\rm \Delta} S_v=5\ {\rm g}\ {\rm kg}^{-1}$, the physical explanation of the observed layer formation is as follows: after the ice starts to melt, there is an accumulation of cold water outside the thin fresh boundary layer. Being heavier than the surrounding fluid, the cold water sinks. The surrounding fluid, however, becomes denser with depth as the local salinity also increases, and eventually the cold water reaches a neutral buoyancy level, thus producing a sequence of layers as seen in figure 3(*a*). In this case, vertically stacked convection rolls form. The layered convection rolls sculpt a layered pattern in the melt front, since these rolls bring warm water from the bulk to the melt front and cause non-uniform heat flux at the interface. At ${\rm \Delta} S_v=10\ {\rm g}\ {\rm kg}^{-1}$, more layers in flow structure appear as compared to the case of ${\rm \Delta} S_v=5\ {\rm g}\ {\rm kg}^{-1}$, in accordance with our physical explanation of neutral buoyancy levels, while the amplitude of the layers becomes smaller than in the case of ${\rm \Delta} S_v=5\ {\rm g}\ {\rm kg}^{-1}$, owing to the weaker flow by the stable stratification. In contrast, at the top, the flow is still stronger because the water is relatively fresh. Thus, the ice morphology is closely linked to flow structures, such as layer formation and convective circulation. The ice morphology regimes depend intrinsically on the flow regimes. Chong *et al.* (Reference Chong, Yang, Wang, Verzicco and Lohse2020) investigated in more detail the flow structures and regime transitions between the convection-dominated regime and the stratification-dominated regime.

We next quantitatively describe the layer height $h$ of the melt front, which occurs in the presence of a stable stratification. We measure the layer thicknesses by calculating the distances between the maximum points at the contour of the ice front ($\phi =0.5$) and then average all these distances as the layer thicknesses. From the analysis of Huppert & Turner (Reference Huppert and Turner1980), by balancing the horizontal density difference due to temperature and the vertical density gradient due to salinity, one can quantify the thickness of these layers as

where $\rho (T,S)$ is the fluid density at temperature $T$ and salinity $S$, ${\rm d}\rho /{\rm d}z$ is the ambient density stratification, and the subscript $\infty$ in $T_\infty$ and $S_\infty$ relates to the mean far-field value.

In figure 3(*b*), we plot the mean layer thickness $h$, normalized by the horizontal density difference, as a function of the stratification. We present data both from our simulations and also from the experiments of Huppert & Turner (Reference Huppert and Turner1980). Our simulation results agree very well quantitatively with the experimental data, with (3.1) also giving a good prediction of the layer thickness $h$. Moreover, our results on the melt front shape match well with those of Huppert & Turner (Reference Huppert and Turner1980), which can also be regarded as validation of our simulations of ice melting in saline water.

## 4. Physical mechanism of the dependence of the average melt rate on salinity

Our objective now is to quantify how salinity variations affect the average melt rate. In figure 4(*a*), we plot the volume of ice $V(t)$ (normalized by its initial value $V_0$) as a function of time for different $S_m$ with ${\rm \Delta} S_v=5\ {\rm g}\ {\rm kg}^{-1}$, ${\rm \Delta} T=20\ {\rm K}$ and $H=5\ {\rm cm}$. Interestingly, the melt rate shows a non-monotonic relationship with the mean ambient salinity $S_m$; see figure 4(*b*). In that figure, to further quantify the melt rate, we have calculated the average melt rate $\bar {f}=1/t_{1/2}$, where $t_{1/2}$ represents the time needed to melt half of the initial volume, shown as a dashed line in figure 4(*a*), and we show $\bar {f}$ as a function of $S_m$ for various ${\rm \Delta} S_v$. We normalize these results relative to the melt rate observed in the absence of salinity. From the data points, one can see a non-monotonic dependence of $\bar {f}$ on $S_m$ in both 2D and 3D simulations: as $S_m$ increases, $\bar {f}$ first decreases and then increases, with a local minimum point depending on ${\rm \Delta} S_v$. Note that we also tried using different thresholds to calculate $\bar {f}$; these change the absolute value of $\bar {f}$, but the trend remains the same. The drivers of this non-monotonicity are not obvious, motivating us to analyse what causes a minimum in the melt rate.

To understand this behaviour, we consider the force (per unit mass) balance in the system. A similar argument was adopted for melting in fresh water (Yang *et al.* Reference Yang, Chong, Liu, Verzicco and Lohse2022*a*), where the density maximum plays an important role. An illustration of the main sources of dominating terms driving the flow is shown in figure 5(*a*). The work of raising/sinking the fluid parcel in the stable stratification (${\rm \Delta} S_v$) is done by the buoyancy force, which is driven by both temperature ${\rm \Delta} T$ and salinity $S_m$. When $S_m$ is small, the temperature dominates the buoyancy, and the cold fresh meltwater moves downward (since the temperature is much larger than $4\,^\circ$C, so that the density maximum effect can be neglected). A circulation flow forms and transports cold water away from the ice and warm water towards the ice. When $S_m$ is large, the salinity dominates the buoyancy, and the cold fresh meltwater moves upward; a circulation flow is again generated and efficiently melts the ice. However, at intermediate $S_m$, the temperature- and salinity-induced buoyancy compensate for each other. In this case, fresh meltwater has almost the same density as the surroundings, which weakens the buoyancy-driven flow. The weakened flow near the melt front can be seen from figure 5(*b*), where we plot the vertical velocity profile in the horizontal direction for three different values of $S_m$. The value $S_m = 3.5\ {\rm g}\ {\rm kg}^{-1}$ is the one where the buoyancy due to the salt concentration and the buoyancy due to thermal effects compensate for each other, and where hence we find the lowest melting rate. More details of the velocity fields are shown in figure 6, following the same cases shown in figure 5(*b*). Frome figure 6(*b*,*c*), one can see that the magnitude of the velocity is correlated with the local melt rates: the location with a faster melt rate also has a larger velocity.

Quantitatively, when the stable stratification is weak, we consider the force balance between the thermal-driven buoyancy to raise the fluid parcel $F_{{\rm \Delta} T}=\frac {1}{2}C_bg{\rm \Delta} T^2$, and the saline-driven buoyancy $F_{{S_m}}=b_0gS_m$, both based on the equation of state (2.8) given above. We obtain

which means that the temperature- and salinity-induced buoyancies compensate for each other.

To check whether the results agree with the theory, we calculate the values of $\varLambda _T$ and $\varLambda _S$ corresponding to all minimum melt rate points. Then, in figure 5(*c*), we plot them in the parameter space spanned by $\varLambda _T$ and $\varLambda _S$ for various ${\rm \Delta} T$ and $H$. It can be seen that the data points for different ${\rm \Delta} T$ and $H$ follow the same trend of $\varLambda _T=1$ (I–II) at low ${\rm \Delta} S_v$. We also notice that at high ${\rm \Delta} S_v$, the minimum melt rate points instead follow roughly $\varLambda _S\sim 1$ (II–III) and are independent of the $\varLambda _T$, with a transition between the different regimes. We give our explanation on this second regime (II–III): when the stable stratification is strong (e.g. for large ${\rm \Delta} S_v$, rightmost column of figure 3*a*), at low $S_m$ a layer of fresh water emerges at the top, which has melted away faster because of the reduced local stratification. At large $S_m$, the salinity-induced buoyancy becomes larger, also resulting in stronger convection and thus faster melting. Therefore, the melt rate has a minimum at medium $S_m$, which results from the competition between buoyancy-driving and local stratification, which is modified by the fresh layer formation at low $S_m$.

## 5. Conclusions and outlook

In summary, we have numerically studied ice melting in saline water, using direct numerical simulations with a realistic, nonlinear equation of state. We have shown that the dependence of the melt rate on the ambient salinity is non-monotonic: as the ambient mean salinity increases, the melt rate first decreases and then increases. The physical origin of this non-monotonic dependence is the competition between thermally-driven buoyancy and salt-driven buoyancy, as well as the stable stratification due to the vertical salinity gradient. We have derived a theoretical model based on a force balance, which collapses the points of the minimum melt rate in the non-monotonic trend. Finally, we have shown the effect of salinity on the melt shape. Layered structures on the melt front have been observed, and a comparison of the layer thicknesses with the experimental results of Huppert & Turner (Reference Huppert and Turner1980) shows quantitative agreement.

From a broader perspective, our results show that the phase-field method is capable of quantitatively modelling the melting process in multi-component turbulent flows (Hester *et al.* Reference Hester, Couston, Favier, Burns and Vasil2020). Our numerical and theoretical results on the effect of salinity on the ice melt rate and shape can be applied to various saline water systems of geophysical relevance, e.g. sea ice, ice shelves, and icy moons (which usually have even higher salinity than we have explored), and more generally to multi-component phase change materials. Future work may address in more detail the effect of salinity on melting point and may enable closer comparison to geophysical and astrophysical ocean environments where lower temperatures, higher salinities and larger scales are relevant. Because large salinities require extremely high resolution, it is difficult to simulate a realistic open ocean environment (${\sim }35\ {\rm g}\ {\rm kg}^{-1}$) with direct numerical simulations. At these salinity values, the density maximum in temperature would disappear (see figure 2*b*), with salinity differences being the main driver of buoyancy. In cases where low salinity remains relevant, a lower ambient temperature would also amplify the density maximum effect. The thermally-driven circulation may change appreciably in this regime within layers of near-uniform salinity.

## Funding

We acknowledge PRACE for awarding us access to MareNostrum in Spain at the Barcelona Computing Center (BSC) under the projects 2020235589 and 2021250115, as well as the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC). We also acknowledge support from the Deutsche Forschungsgemeinschaft Priority Programme SPP 1881, ‘Turbulent Superstructures’. This research was supported in part by the National Science Foundation under grant no. NSF PHY-1748958.

## Declaration of interests

The authors report no conflict of interest.