Hostname: page-component-5d59c44645-mrcq8 Total loading time: 0 Render date: 2024-02-24T10:04:25.451Z Has data issue: false hasContentIssue false

Ice melting in salty water: layering and non-monotonic dependence on the mean salinity

Published online by Cambridge University Press:  22 August 2023

Rui Yang*
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands
Christopher J. Howland*
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands
Hao-Ran Liu*
Affiliation:
Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, PR China
Roberto Verzicco*
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands Dipartimento di Ingegneria Industriale, University of Rome ‘Tor Vergata’, Roma 00133, Italy Gran Sasso Science Institute - Viale F. Crispi, 7 67100 L'Aquila, Italy
Detlef Lohse*
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands Max Planck Institute for Dynamics and Self-Organization, 37077 Göttingen, Germany

Abstract

The presence of salt in seawater strongly affects the melt rate and the shape evolution of ice, both of utmost relevance in ice–ocean interactions and thus for the climate. To get a better quantitative understanding of the physical mechanics at play in ice melting in salty water, we numerically investigate the lateral melting of an ice block in stably stratified saline water. The developing ice shape from our numerical results shows good agreement with the experiments and theory from Huppert & Turner (J. Fluid Mech., vol. 100, 1980, pp. 367–384). Furthermore, we find that the melt rate of ice depends non-monotonically on the mean ambient salinity: it first decreases for increasing salt concentration until a local minimum is attained, and then increases again. This non-monotonic behaviour of the ice melt rate is due to the competition among salinity-driven buoyancy, temperature-driven buoyancy and salinity-induced stratification. We develop a theoretical model based on the force balance which gives a prediction of the salt concentration for which the melt rate is minimal, and is consistent with our data. Our findings give insight into the interplay between phase transitions and double-diffusive convective flows.

Type
JFM Rapids
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press.

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 Sun2021a) and Wang et al. (Reference Wang, Calzavarini, Sun and Toschi2021b,Reference Wang, Jiang, Du, Sun and Calzavarinic) 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 Lohse2022a). 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

(2.1)\begin{gather} \partial_t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\rho_0^{{-}1} \boldsymbol{\nabla} p + \nu \nabla^2 \boldsymbol{u} - (g\rho'/\rho_0) \boldsymbol{e_z} - \phi \boldsymbol{u} / \eta , \end{gather}
(2.2)\begin{gather}\partial_t T + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \kappa_T \nabla^2 T + (\mathcal{L}/c_p) \partial_t \phi , \end{gather}
(2.3)\begin{gather}\partial_t S + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} S = \kappa_S \nabla^2 S + (S\partial_t \phi - \kappa_S\boldsymbol{\nabla} \phi \boldsymbol{\cdot} \boldsymbol{\nabla} S)/(1 - \phi + \delta), \end{gather}
(2.4)\begin{gather}\partial_t \phi = D \nabla^2 \phi - \frac{D}{\varepsilon^2} \phi(1-\phi)\left(1-2\phi + \frac{\varepsilon}{\gamma}(T + m S)\right) . \end{gather}

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$:

(2.5)\begin{equation} \mathcal{L} v_n = c_p\kappa_T ( \partial_n T^{(S)}-\partial_n T^{(L)}), \kappa_S\partial_n S^{(L)}+S^{(L)}v_n=0. \end{equation}

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:

(2.6a,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:

(2.7ac)\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

(2.8)\begin{equation} \rho'={-C_b}/{2}(T-T_0-c_S S)^2+b_0 S,\end{equation}

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$:

(2.9ac)\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}
(2.10ac)\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$:

(2.11a,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 Caulfield2022b). 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 Lohse2022a).

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.

Figure 1. (a) The normalized volume of ice as a function of time for the case of $S_m=5\ {\rm g}\ {\rm kg}^{-1}$ and ${\rm \Delta} S_v=5\ {\rm g}\ {\rm kg}^{-1}$. Different base resolutions are tested (the refined resolution is fixed as five times larger than the base resolution). (b) The overall melt rate (with the same definition as in the main text) as a function of base resolution. Convergence is shown as the resolution increases. In this case, we choose $n_x=288$.

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(ce) 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$.

Figure 2. (a) An illustration of the simulation set-up. Initially, the ice block is set at the right sidewall, the temperature of the water is set to be uniform as ${\rm \Delta} T$ and the salinity of the water is set with a vertical gradient depending on $S_m$ and ${\rm \Delta} S_v$. (b) Plot of relative density $\varDelta \rho$ versus temperature for different salinities. The dots mark the melting points. (ce) Snapshots of temperature (left column), salinity field (centre) and contour of melt front (right) at (c) $S_v=0$, (d) $S_v=5\ {\rm g}\ {\rm kg}^{-1}$ and (e) $S_v=10\ {\rm g}\ {\rm kg}^{-1}$.

At relatively low or zero vertical salinity variation ${\rm \Delta} S_v =0$ (figure 2c), 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 2d), 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 2e), 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.

Figure 3. (a) Instantaneous snapshots of temperature field (left) and salinity field (right) for various $S_m$ and ${\rm \Delta} S_v$. The velocity field is shown by arrows in the temperature field. (b) The layer thickness $h$ per unit density as function of the density gradient, following the same representation as in figure 10 of Huppert & Turner (Reference Huppert and Turner1980). The reference density is $\rho _0$ and $\varDelta \rho _T=\rho (0, S_{\infty })-\rho (T_{\infty }, S_{\infty })$. The dashed line represents the theoretical result ($h$) (3.1).

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

(3.1)\begin{equation} h=(0.65 \pm 0.06)[\rho(0, S_{\infty})-\rho(T_{\infty}, S_{\infty})]\left(\frac{\mathrm{d} \rho}{\mathrm{d} z}\right)^{{-}1},\end{equation}

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.

Figure 4. (a) The normalized volume of ice $V(t)/V_0$ as a function of time $t/t_f$, where $t_f$ is the free-fall time unit. The dashed line represents the location of half of the initial ice volume. (b) The normalized melt rate $\bar {f}/\bar {f}_0$ as a function of $S_m$ for various ${\rm \Delta} S_v$ (colour-coded), where $\bar {f}_0$ is the melt rate without salinity. Circle data points represent 2D simulations, and square data points represent 3D simulations. The red circle data points represent the location of the minimum $\bar {f}$.

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

Figure 5. (a) An illustration of the effects of temperature and salinity. The red and blue colours represent the buoyancy force driven by $T$ and $S$, with the arrows showing the direction of buoyancy. The black line represents the stable stratification. (b) The instantaneous vertical velocity profile (in free-fall velocity units) at mid-height along $x$. The quantity $S_m=3.5\ {\rm g}\ {\rm kg}^{-1}$ corresponds to the minimum melt rate. The inset figure illustrates the $T(x)$, $S(x)$ and $\rho (x)$ profiles at different ambient $S_m$. The dashed line shows the location of the melt front. (c) Location of the minimal melt rate in the parameter space spanned by $\varLambda _T$ and $\varLambda _S$ for various ${\rm \Delta} T$ and $H$. The dashed lines show the prediction from (4.1) and $\varLambda _S=1$ for reference. Regime I is ‘$T$-driven buoyancy’, regime II is ‘$S$-driven buoyancy’ and regime III is ‘fresh layer formation’. The shaded area shows the upper and lower bounds of $\varLambda _S$ from the simulation results.

Figure 6. The instantaneous vertical velocity profiles at different heights along $x$. for ${\rm \Delta} S_v=0$ and (a) $S_m=0$; (b) $S_m=3.5\ {\rm g}\ {\rm kg}^{-1}$; (c) $S_m=10\ {\rm g}\ {\rm kg}^{-1}$, the same cases as in figure 5(b). (df) shows the instantaneous vertical velocity profiles at mid-height as a function of time, corresponding to the cases in (ac). The black lines show the location of the ice front.

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

(4.1)\begin{equation} \tfrac{1}{2}C_bg{\rm \Delta} T^2=b_0gS_m\quad {\rm or}\quad \varLambda_T=1,\end{equation}

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 3a), 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 2b), 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.

References

Cao, Y., Tan, W. & Wu, Z. 2018 Aircraft icing: an ongoing threat to aviation safety. Aerosp. Sci. Technol. 75, 353385.CrossRefGoogle Scholar
Cenedese, C. & Straneo, F. 2023 Icebergs melting. Annu. Rev. Fluid Mech. 55, 377402.CrossRefGoogle Scholar
Chong, K.L., Yang, R., Wang, Q., Verzicco, R. & Lohse, D. 2020 Café latte: spontaneous layer formation in laterally cooled double diffusive convection. J. Fluid Mech. 900, R6.CrossRefGoogle Scholar
Couston, L.-A., Hester, E., Favier, B., Taylor, J.R., Holland, P.R. & Jenkins, A. 2021 Topography generation by melting and freezing in a turbulent shear flow. J. Fluid Mech. 911, A44.CrossRefGoogle Scholar
Davis, S.H., Müller, U. & Dietsche, C. 1984 Pattern selection in single-component systems coupling Bénard convection and solidification. J. Fluid Mech. 144, 133151.CrossRefGoogle Scholar
Dhaidan, N.S. & Khodadadi, J. 2015 Melting and convection of phase change materials in different shape containers: a review. Renew. Sustain. Energy Rev. 43, 449477.CrossRefGoogle Scholar
Du, Y., Wang, Z., Jiang, L., Calzavarini, E. & Sun, C. 2023 Sea water freezing modes in a natural convection system. J. Fluid Mech. 960, A35.CrossRefGoogle Scholar
Edwards, T.L., et al. 2021 Projected land ice contributions to twenty-first-century sea level rise. Nature 593 (7857), 7482.CrossRefGoogle ScholarPubMed
Favier, B., Purseed, J. & Duchemin, L. 2019 Rayleigh–Bénard convection with a melting boundary. J. Fluid Mech. 858, 437473.CrossRefGoogle Scholar
Hester, E.W., Couston, L.-A., Favier, B., Burns, K.J. & Vasil, G.M. 2020 Improved phase-field models of melting and dissolution in multi-component flows. Proc. R. Soc. A 476 (2242), 20200508.CrossRefGoogle ScholarPubMed
Hester, E.W., McConnochie, C.D., Cenedese, C., Couston, L.-A. & Vasil, G. 2021 Aspect ratio affects iceberg melting. Phys. Rev. Fluids 6 (2), 023802.CrossRefGoogle Scholar
Holland, M.M., Bitz, C.M. & Tremblay, B. 2006 Future abrupt reductions in the summer arctic sea ice. Geophys. Res. Lett. 33 (23), L23503.CrossRefGoogle Scholar
Howland, C.J. 2022 AFiD-MuRPhFi Documentation. Available at: https://chowland.github.io/AFiD-MuRPhFi/.Google Scholar
Huppert, H.E. & Turner, J.S. 1978 On melting icebergs. Nature 271 (5640), 4648.CrossRefGoogle Scholar
Huppert, H.E. & Turner, J.S. 1980 Ice blocks melting into a salinity gradient. J. Fluid Mech. 100 (2), 367384.CrossRefGoogle Scholar
Kang, W. & Flierl, G. 2020 Spontaneous formation of geysers at only one pole on Enceladus's ice shell. Proc. Natl Acad. Sci. 117 (26), 1476414768.CrossRefGoogle ScholarPubMed
Kang, W., Mittal, T., Bire, S., Campin, J.-M. & Marshall, J. 2022 How does salinity shape ocean circulation and ice geometry on enceladus and other icy satellites? Sci. Adv. 8 (29), eabm4665.CrossRefGoogle ScholarPubMed
Kooij, G.L., Botchev, M.A., Frederix, E.M., Geurts, B.J., Horn, S., Lohse, D., van der Poel, E.P., Shishkina, O., Stevens, R.J.A.M. & Verzicco, R. 2018 Comparison of computational codes for direct numerical simulations of turbulent Rayleigh–Bénard convection. Comput. Fluids 166, 18.CrossRefGoogle Scholar
Kranenborg, E.J. & Dijkstra, H.A. 1998 On the evolution of double-diffusive intrusions into a stably stratified liquid: a study of the layer merging process. Intl J. Heat Mass Transfer 41 (18), 27432756.CrossRefGoogle Scholar
Liu, H.-R., Ng, C.S., Chong, K.L., Lohse, D. & Verzicco, R. 2021 An efficient phase-field method for turbulent multiphase flows. J. Comput. Phys. 446, 110659.CrossRefGoogle Scholar
McConnochie, C.D. & Kerr, R.C. 2016 The effect of a salinity gradient on the dissolution of a vertical ice face. J. Fluid Mech. 791, 589607.CrossRefGoogle Scholar
Ostilla-Mónico, R., Yang, Y., Van Der Poel, E.P., Lohse, D. & Verzicco, R. 2015 A multiple-resolution strategy for direct numerical simulation of scalar turbulence. J. Comput. Phys. 301, 308321.CrossRefGoogle Scholar
van der Poel, E.P., Ostilla-Mónico, R., Donners, J. & Verzicco, R. 2015 A pencil distributed finite difference code for strongly turbulent wall-bounded flows. Comput. Fluids 116, 1016.CrossRefGoogle Scholar
Purseed, J., Favier, B., Duchemin, L. & Hester, E.W. 2020 Bistability in Rayleigh-Bénard convection with a melting boundary. Phys. Rev. Fluids 5 (2), 023501.CrossRefGoogle Scholar
Ravichandran, S., Toppaladoddi, S. & Wettlaufer, J.S. 2022 The combined effects of buoyancy, rotation, and shear on phase boundary evolution. J. Fluid Mech. 941, A39.CrossRefGoogle Scholar
Ravichandran, S. & Wettlaufer, J.S. 2021 Melting driven by rotating Rayleigh–Bénard convection. J. Fluid Mech. 916, A28.CrossRefGoogle Scholar
Ristroph, L. 2018 Sculpting with flow. J. Fluid Mech. 838, 14.CrossRefGoogle Scholar
Roquet, F., Madec, G., Brodeau, L. & Nycander, J. 2015 Defining a simplified yet ‘realistic’ equation of state for seawater. J. Phys. Oceanogr. 45 (10), 25642579.CrossRefGoogle Scholar
Rubinstein, L. 1971 The Stefan Problem, Translations of Mathematical Monographs, vol. 27. American Mathematical Society.Google Scholar
Spencer, J., Pearl, J., Segura, M., Flasar, F., Mamoutkine, A., Romani, P., Buratti, B., Hendrix, A., Spilker, L. & Lopes, R. 2006 Cassini encounters enceladus: background and the discovery of a south polar hot spot. Science 311 (5766), 14011405.CrossRefGoogle ScholarPubMed
Toppaladoddi, S. 2021 Nonlinear interactions between an unstably stratified shear flow and a phase boundary. J. Fluid Mech. 919, A28.CrossRefGoogle Scholar
Wang, Z., Calzavarini, E. & Sun, C. 2021 a Equilibrium states of the ice-water front in a differentially heated rectangular cell (a). Europhys. Lett. 135 (5), 54001.CrossRefGoogle Scholar
Wang, Z., Calzavarini, E., Sun, C. & Toschi, F. 2021 b How the growth of ice depends on the fluid dynamics underneath. Proc. Natl Acad. Sci. 118, 10.Google ScholarPubMed
Wang, Z., Jiang, L., Du, Y., Sun, C. & Calzavarini, E. 2021 c Ice front shaping by upward convective current. Phys. Rev. Fluids 6 (9), L091501.CrossRefGoogle Scholar
Wang, S.L., Sekerka, R.F., Wheeler, A.A., Murray, B.T., Coriell, S.R., Braun, R.J. & McFadden, G.B. 1993 Thermodynamically-consistent phase-field models for solidification. Physica D 69 (1–2), 189200.CrossRefGoogle Scholar
Weady, S., Tong, J., Zidovska, A. & Ristroph, L. 2022 Anomalous convective flows carve pinnacles and scallops in melting ice. Phys. Rev. Lett. 128 (4), 044502.CrossRefGoogle ScholarPubMed
Worster, M.G. 1997 Convection in mushy layers. Annu. Rev. Fluid Mech. 29 (1), 91122.CrossRefGoogle Scholar
Yang, R., Chong, K.L., Liu, H.-R., Verzicco, R. & Lohse, D. 2022 a Abrupt transition from slow to fast melting of ice. Phys. Rev. Fluids 7 (8), 083503.CrossRefGoogle Scholar
Yang, R., Howland, C.J., Liu, H.-R., Verzicco, R. & Lohse, D. 2023 Morphology evolution of a melting solid layer above its melt heated from below. J. Fluid Mech. 956, A23.CrossRefGoogle Scholar
Yang, Y., Verzicco, R. & Lohse, D. 2016 From convection rolls to finger convection in double-diffusive turbulence. Proc. Natl Acad. Sci. 113, 6973.CrossRefGoogle ScholarPubMed
Yang, Y., Verzicco, R., Lohse, D. & Caulfield, C. 2022 b Layering and vertical transport in sheared double-diffusive convection in the diffusive regime. J. Fluid Mech. 933, A30.CrossRefGoogle Scholar
Figure 0

Figure 1. (a) The normalized volume of ice as a function of time for the case of $S_m=5\ {\rm g}\ {\rm kg}^{-1}$ and ${\rm \Delta} S_v=5\ {\rm g}\ {\rm kg}^{-1}$. Different base resolutions are tested (the refined resolution is fixed as five times larger than the base resolution). (b) The overall melt rate (with the same definition as in the main text) as a function of base resolution. Convergence is shown as the resolution increases. In this case, we choose $n_x=288$.

Figure 1

Figure 2. (a) An illustration of the simulation set-up. Initially, the ice block is set at the right sidewall, the temperature of the water is set to be uniform as ${\rm \Delta} T$ and the salinity of the water is set with a vertical gradient depending on $S_m$ and ${\rm \Delta} S_v$. (b) Plot of relative density $\varDelta \rho$ versus temperature for different salinities. The dots mark the melting points. (ce) Snapshots of temperature (left column), salinity field (centre) and contour of melt front (right) at (c) $S_v=0$, (d) $S_v=5\ {\rm g}\ {\rm kg}^{-1}$ and (e) $S_v=10\ {\rm g}\ {\rm kg}^{-1}$.

Figure 2

Figure 3. (a) Instantaneous snapshots of temperature field (left) and salinity field (right) for various $S_m$ and ${\rm \Delta} S_v$. The velocity field is shown by arrows in the temperature field. (b) The layer thickness $h$ per unit density as function of the density gradient, following the same representation as in figure 10 of Huppert & Turner (1980). The reference density is $\rho _0$ and $\varDelta \rho _T=\rho (0, S_{\infty })-\rho (T_{\infty }, S_{\infty })$. The dashed line represents the theoretical result ($h$) (3.1).

Figure 3

Figure 4. (a) The normalized volume of ice $V(t)/V_0$ as a function of time $t/t_f$, where $t_f$ is the free-fall time unit. The dashed line represents the location of half of the initial ice volume. (b) The normalized melt rate $\bar {f}/\bar {f}_0$ as a function of $S_m$ for various ${\rm \Delta} S_v$ (colour-coded), where $\bar {f}_0$ is the melt rate without salinity. Circle data points represent 2D simulations, and square data points represent 3D simulations. The red circle data points represent the location of the minimum $\bar {f}$.

Figure 4

Figure 5. (a) An illustration of the effects of temperature and salinity. The red and blue colours represent the buoyancy force driven by $T$ and $S$, with the arrows showing the direction of buoyancy. The black line represents the stable stratification. (b) The instantaneous vertical velocity profile (in free-fall velocity units) at mid-height along $x$. The quantity $S_m=3.5\ {\rm g}\ {\rm kg}^{-1}$ corresponds to the minimum melt rate. The inset figure illustrates the $T(x)$, $S(x)$ and $\rho (x)$ profiles at different ambient $S_m$. The dashed line shows the location of the melt front. (c) Location of the minimal melt rate in the parameter space spanned by $\varLambda _T$ and $\varLambda _S$ for various ${\rm \Delta} T$ and $H$. The dashed lines show the prediction from (4.1) and $\varLambda _S=1$ for reference. Regime I is ‘$T$-driven buoyancy’, regime II is ‘$S$-driven buoyancy’ and regime III is ‘fresh layer formation’. The shaded area shows the upper and lower bounds of $\varLambda _S$ from the simulation results.

Figure 5

Figure 6. The instantaneous vertical velocity profiles at different heights along $x$. for ${\rm \Delta} S_v=0$ and (a) $S_m=0$; (b) $S_m=3.5\ {\rm g}\ {\rm kg}^{-1}$; (c) $S_m=10\ {\rm g}\ {\rm kg}^{-1}$, the same cases as in figure 5(b). (df) shows the instantaneous vertical velocity profiles at mid-height as a function of time, corresponding to the cases in (ac). The black lines show the location of the ice front.