Sea water freezing modes in a natural convection system

Abstract

Sea ice is crucial in many natural processes and human activities. Understanding the dynamical couplings between the inception, growth and equilibrium of sea ice and the rich fluid mechanical processes occurring at its interface and interior is of relevance in many domains, ranging from geophysics to marine engineering. Here we investigate experimentally the complete freezing process of water with dissolved salt in a standard natural convection system, i.e. the prototypical Rayleigh–Bénard cell. Due to the presence of a mushy phase, the studied system is considerably more complex than the freezing of fresh water in the same conditions (Wang et al., Proc. Natl Acad. Sci. USA, vol. 118, issue 10, 2021, e2012870118). We measure the ice thickness and porosity at the dynamical equilibrium state for different initial salinities of the solution and temperature gaps across the cell. These observables are non-trivially related to the controlling parameters of the system as they depend on the heat transport mode across the cell. We identify in the experiments five out of the six possible modes of heat transport. We highlight the occurrence of brine convection through the mushy ice and of penetrative convection in stably stratified liquid underlying the ice. A one-dimensional multi-layer heat flux model built on the known scaling relations of global heat transport in natural convection systems in liquids and porous media is proposed. Given the measured porosity of the ice, it allows us to predict the corresponding ice thickness, in a unified framework.

1. Introduction

The evolution of sea ice has important impacts on many environmental and geological processes as well as human activities. Examples include ocean circulation (Clark, Alley & Pollard 1999; Joughin, Alley & Holland 2012; Hanna et al. 2013; Straneo & Heimbach 2013; Stevens et al. 2020), global sea-level rise (Wadhams & Munk 2004), land-surface albedo (Curry, Schramm & Ebert 1995; Perovich et al. 2002; Scagliarini et al. 2020), biodiversity (Post et al. 2013), microplastic dispersion and sequestration (Peeken et al. 2018; Obbard et al. 2014), and winter navigation at high latitudes and polar areas (de Andrés, Saarinen & Uuskallio 2018; Ho 2010). Generally, sea ice evolution is associated with the complex fluid dynamic processes occurring in the oceans that involve wide ranges of spatial and temporal scales. Their descriptions through numerical models in the geophysical context make use of so-called micro-scale physical parametrizations that have a key role in the reliability and accuracy of the resulting predictions (Rae et al. 2015; Rousset et al. 2015).

In this study, we focus on the role of buoyancy-driven natural convection in determining upper-surface sea ice growth and the attainment of dynamical equilibrium. We do this by performing experiments at the laboratory scale in a highly controllable setting, the Rayleigh–Bénard cell (see below). With this approach, we aim at providing a highly reliable small-scale parametrization of the process that is crucial for geophysical-scale computational models in the wide range of applications mentioned above.

The Rayleigh–Bénard (RB) system, a fluid layer parallelly confined between a warmer horizontal bottom plate and a cooler top plate, and insulated on the lateral sides (Ahlers, Grossmann & Lohse 2009; Lohse & Xia 2010; Chillà & Schumacher 2012), is a paradigmatic model system in the physics of nonlinear systems and fluid mechanics. Its behaviour has been studied for decades, and it has been crucial in developing the understanding of hydrodynamic instabilities, transition to chaos, turbulence and turbulent transport. Its local and global properties across its wide spectrum of regimes are nowadays well-characterized (Ahlers et al. 2009; Chillà & Schumacher 2012). More recently, the RB system has been proven to be also convenient for studying the couplings between liquid–solid phase change and fluid flow phenomena (Esfahani et al. 2018; Favier, Purseed & Duchemin 2019). It has a laboratory scale, it is closed, and in/out energy fluxes can be measured, hence it is highly controllable, yet able to reproduce the rich process of water freezing in the presence of unsteady and turbulent flow conditions. Several recent studies have investigated the influence of temperature (Wang et al. 2021c; Yang et al. 2023b), inclination (Wang et al. 2021d; Yang et al. 2022), rotation (Ravichandran, Toppaladoddi & Wettlaufer 2022), initial conditions, and aspect ratio (Wang, Calzavarini & Sun 2021b) on fresh water freezing/melting in the RB system. However, the presence of dissolved salt in sea water adds to the complexity of the processes. In particular, two factors need to be considered to properly study the evolution of sea ice.

The first factor is the mushy structure of sea ice (Feltham et al. 2006), which allows for the occurrence of interstitial fluid motion, i.e. brine drainage and thermal convection across the sea ice. Various studies have discussed the mushy structure forming in the solidification of aqueous solutions. The focus has mainly been on onset conditions for convective brine drainage (Wettlaufer, Worster & Huppert 1997), mushy phase morphology (Peppin, Huppert & Worster 2008), ice layer growth rate (Kerr et al. 1990), ice porosity (Wettlaufer et al. 1997), flow structure within the mushy phase (Chen 1995; Worster & Wettlaufer 1997), and salinity in the liquid phase (Wettlaufer et al. 1997; Peppin et al. 2008). While these studies together with relevant reviews (Worster 1997; Wells, Hitchen & Parkinson 2019; Anderson & Guba 2020) provide a detailed, thorough and multi-perspective understanding of brine drainage and mushy phase structure, the interactions between the mushy structure and a neighbouring vigorous thermal convective flow are much less explored. A few studies discussing the growth of sea ice focus either on establishing governing momentum and energy equations suitable for numerical simulations (Feltham et al. 2006; Covello et al. 2016; Wells et al. 2019), or on theoretical models that view ice growth as purely diffusion controlled (Notz 2005; Worster & Rees Jones 2015).

The second factor is salty water's peculiar equation of state, in particular its density dependence on salinity and temperature (Gebhart & Mollendorf 1977), and its freezing point depression (Bodnar 1993). While for normal sea water with 3.5 % salinity (i.e. the mass concentration of dissolved salt), the density decreases monotonically as the temperature is increased from its freezing point, this is not the case for less salty water (encountered, for instance, near river estuaries) or fresh water, where a maximal density occurs in the interval $[0,4]\,^{\circ }$C (Lyman & Fleming 1940; Pawlowicz 2015). Salty water's peculiar density dependence may lead to different combinations of stable/unstable stratifications of the liquid phase in the RB cell. In a similar yet simpler scenario, the combination of stable stratification and unstably stratified turbulence has been found to have a significant influence on the freezing dynamics and equilibrium ice thickness of fresh water (Townsend 1964; Wang et al. 2021c). Other studies point out that the kind of stratification, either stable or unstable, in touch with the melting front has a key role in determining the melting rate of a pure-ice layer in aqueous solutions (Sugawara & Irvine 2000; Mondal et al. 2019; Rosevear, Gayen & Galton-Fenzi 2021; Yang et al. 2023a). However, these studies focus on stratifications due to the salinity gradient in double-diffusive convection with linearized equation of state salt–water mixtures, i.e. neglecting the above-mentioned important nonlinear dependencies of density with respect to temperature and salinity.

In this work, we investigate how ice growth and fluid dynamics in salty water couple and co-evolve until the equilibrium in the RB natural convection system is reached. We aim to reveal how sea ice growth rate depends on controlled environmental conditions through its interactions with mushy phase convection and with the variety of density stratification patterns and fluid flows occurring in the underlying liquid phase.

The paper is organized as follows. Section 2 describes the experimental set-up and the measurement methodology. Section 3 describes the experimental results, focusing on the growth dynamics and equilibrium morphology of the ice layer. To provide a first-step understanding of the complex process, in § 4 we discuss how the equilibrium ice thickness, its porosity and the saturation time depend on temperature and salinity conditions, as well as the correlation among these global properties. Section 5 introduces a one-dimensional multi-layer heat flux model that describes the equilibrium ice thickness's dependence on the system control parameters, and identifies five out of six possible modes of heat transport in the experimental system. Finally, we end with conclusions and outlook in § 6.

2. Experiments

The experiments are performed in a rectangular cuboid RB cell (height $H=0.12$ m, aspect ratio $\varGamma = L/H = 2$, width–height ratio $W/H = 0.5$, i.e. a so-called quasi-two-dimensional system). An expansion vessel is attached to compensate for the volume change during the freezing process and so to monitor the mushy phase growth; see figure 1(a). The top and bottom plate temperatures are measured with thermistors and controlled with circulating baths. Sidewalls are isolated with foam material to approach ideal adiabatic boundary conditions. The cell is initially filled with degassed salty water of initial salinity $S_i$. Here, we use sodium chloride (NaCl), the principal solute of sea water (Lyman & Fleming 1940), as the salt, and $S_i$ is varied from 0 % (fresh water) to 3.5 % (average ocean salinity), with the corresponding freezing point $T_{0i}$ (Hall, Sterner & Bodnar 1988) varying from $0^{\circ }$C to $-2.1\,^{\circ }$C. The (hot) bottom plate temperature $T_b$, which is equal to the initial liquid temperature, ranges from $0.4^{\circ }$C to 12.9 $^{\circ }$C, and the (cold) top plate temperature $T_t$ is set 10 K below the freezing temperature, hence $T_t = T_{0i}-10$ K (from $-10^{\circ }$C to $-12.1\,^{\circ }$C, depending on $S_i$). A dimensionless bottom superheat parameter can then be defined as $\varTheta _i = (T_b-T_{0i})/(T_b-T_t)$. This parameter characterizes the initial ratio of the admissible temperature difference in the liquid phase over the total gap across the cell. It can also evolve slightly during the experiment. This is because the melting temperature depends on the actual salinity, and at the equilibrium state, $S_e$ and so $T_{0e}$ might differ from the initial values (actually $S_e$ can be as high as $7\,\%$, with $T_{0e}\approx -4.2\,^{\circ }$C; see figure 5b). Since the total quantity of solute is constant, one expects an enhancement of salinity in the liquid due to salt ejection from the mush (and so an increase of $\varTheta$).

Figure 1. Sketch of the experimental system and temperature–salinity parameter space. (a) Cartoon of the experimental set-up. The height of the cell is $H=0.12$ m, the aspect ratio is $\varGamma = L/H = 2$, and the width–height ratio is $W/H = 0.5$. An expansion vessel is attached to compensate for the volume change and to monitor the ice/mushy phase growth. The bottom plate temperature is kept at $T_b$, and the top plate temperature $T_t$ is maintained at 10 K below the freezing point $T_{0i}$ at the initial salinity $S_i$. A dimensionless bottom superheat parameter is defined as $\varTheta _i=(T_b-T_{0i})/(T_b-T_t)$. (b) Temperature–salinity parameter space. In this work, we test $\varTheta _i$ values $1/5$ (purple), $1/3$ (red), $3/8$ (green), $3/7$ (pink), $1/2$ (blue) and $3/5$ (brown). The initial salinity $S_i$ varies from 0 % to 3.5 %; $T_b$ (dots) is shown for each case. We also plot $T_{0i}$ (black squares) and $T_t$ (grey triangles) as functions of $S_i$.

In this work, we first test the states with $\varTheta _i$ values $1/3$, $3/8$ and $1/2$ with salinities $S_i \in [0, 3.5]\%$. Cases with $\varTheta _i$ values $1/5$, $3/7$ and $3/5$, and $S_i = 3.5$%, are also performed to test the results in a wider temperature range. A visualization of the temperature–salinity parameter space explored in the experiments is provided in figure 1(b). It reports the temperature values $T_b$, $T_{0i}$ and $T_t$ for different $\varTheta _i$ and $S_i$. In this experimental configuration, the ice/mushy phase grows from the top plate until the system reaches an equilibrium. We assume that this state is attained when the mean mushy phase thickness $h(t)$ (see figure 1a) changes less than 1 % in eight hours, with the maximum experiment duration being three days. The ice or mushy phase thickness and its porosity $\phi (t)$ (i.e. the volume fraction of residual liquid in the ice or mushy phase) can be obtained, respectively, by identifying the ice–liquid interface in photos of the system, and by using the expansion vessel readings together with the assumption of total mass conservation in the system. The evaporation from the expansion vessel is negligible as it causes only a mass loss less than 0.1 % of the total mass of the system. We refer to Appendix A for more details on the measurement protocols.

3. Mushy phase morphology and growth dynamics

Figures 2(a,c) show photos of the experiment at its equilibrium state in two typical cases ($S_i=0$, 3.5 % and $\varTheta _i=3/8$). One can readily appreciate the different thicknesses of the iced domains (thinner in the salty water case) and the different interface shapes between the two cases. However, their dissimilarity is not only volumetric and morphological, but also structural. To illustrate this, we hollowed out a tiny (5 mm) hemispherical cavity on the top horizontal surface of the ice matrix on which we deposited a drop of red dye. Figures 2(b,d) show two lateral view photos of the stain after 5 min to form its deposition. We observe that when $S_i=0$ (figure 2b), the iced phase is transparent and the dye does not spread into it, indicating the dense structure of the matrix. On the contrary, at $S_i=3.5\,\%$ (figure 2d), the iced phase looks opaque and the dye spreads anisotropically, which highlights a complex porous medium filled with liquid (i.e. a mushy phase structure) with internal vertical channels.

Figure 2. Images of the iced layer together with its internal structure, and evolution of the freezing front. (a–d) Photos of ice/mushy phases without (a,b) and with (c,d) porosity. The experimental conditions are (a,b) $S_i=0\,\%$, and (c,d) $S_i=3.5\,\%$, both with $\varTheta _i = 3/8$. In (b,d), the red dye of the same amount and concentration is deposited into a small hemispherical pit carved at the same location on the upper surface of the corresponding solid matrix in (a,c). The ice in (a,b) is transparent, and the dye does not spread. The mushy phase with porosity in (c,d) is relatively opaque, and the dye spreads anisotropically, possibly indicating the presence of channels. (e) Dynamics of the freezing front. The mushy phase thicknesses $h(t)$, non-dimensionalized by $H$, are plotted as functions of time non-dimensionalized by the diffusive time scale $t^* =H^2/\kappa _{ice}=1.22\times 10^4$ s (where $\kappa _{ice}$ is ice thermal diffusivity) for three experimental runs. (f) The same data as in (e) are compensated by the diffusive behaviour $h(t)/\sqrt {t}$ (normalized by the first data to compare different cases).

It is noteworthy that the different iced phase structures correspond to different freezing dynamics. This is well illustrated in figures 2(e,f), where the temporal growth of the freezing front, $h(t)$, is shown for three sets of experimental conditions, corresponding to pure water, intermediate salinity and sea water salinity, all at $\varTheta _i=3/8$. Though all three cases display an initial period of diffusive ice phase growth, i.e. $h \propto \sqrt {t}$, later the growth rate decreases for the case $S_i=0$, 2 %; while when $S_i=3.5\,\%$, the instantaneous growth rate increases and even exceeds the diffusive growth rate before reaching the equilibrium. It is remarkable also to note that the equilibrium height is not simply proportional to the salinity. The above highlighted temporal evolution of the ice growth in a salty solution depends on the internal structure of the mushy phase that in turn might also change in time (e.g. due to the formation of brine channels Wettlaufer et al. 1997; Worster & Wettlaufer 1997). A quantitative description of this process is challenging because it requires detailed knowledge of the mushy microstructure. The same is true for any attempt to characterize the morphology of the solid–fluid interface. For this reason, in the following we focus on a simpler question: we try to tackle the globally averaged equilibrium state that is reached asymptotically in time. In other words, we aim at understanding the factors that control the average thickness, denoted $h_e$, of the iced mushy phase.

4. Global properties of the equilibrium state

We display in figure 3 two important spatially averaged properties of the mushy phase at the equilibrium (from now on denoted with a subscript $e$): they are the dimensionless thickness $h_e/H$ and the porosity $\phi _e$. Besides this, we report the saturation time (or time to equilibrium) denoted as $t_{e}^{90\,\%}$, here for practical reasons defined as the time it takes to reach 90 % of the equilibrium thickness non-dimensionalized by the diffusive time scale $t^*=H^2/\kappa _{ice}=1.22\times 10^4$ s. The dependence of the thickness $h_e$ as a function of the initial salinity $S_i$ is very different for distinct values of the control parameter $\varTheta _i$ (see figure 3a): $h_e$ exhibits a similar non-monotonic (decrease–increase–decrease) pattern for $\varTheta _i = 1/3, 3/8$, and is monotonically decreasing for $\varTheta _i = 1/2$. These trends are in part reflected on the $\phi _e(S_i)$ dependence, which presents considerable fluctuations for intermediate salinity values, although its overall trend is increasing for all $\varTheta _i$ cases (see figure 3b).

Figure 3. Global properties of the equilibrium state from experiments. The initial superheats are respectively $\varTheta _i$ values $1/3$ (red), $3/8$ (green) and $1/2$ (blue), and the initial salinity $S_i$ varies from 0 % to $3.5\,\%$. The error bars are measured by repeated experiments in two cases with $S_i = 3.5\,\%$ and $\varTheta _i = 1/3$, $1/2$, respectively. (a) The mean mushy phase thickness $h_e$ initially decreases with $S_i$ for all three superheats, and is followed by a sharp increase for the $\varTheta _i = 1/3$ and $3/8$ cases. The shaded areas indicate the minimal and maximal spatial variations of mushy phase thicknesses. The thicknesses are non-dimensionalized by $H$. (b) The mean mushy phase porosity $\phi _e$ exhibits a non-monotonic yet overall increasing trend with $S_i$. (c) Thickness $h_e$ does not exhibit a statistically strong correlation with $\phi _e$. (d) The saturation time $t_e^{90\,\%}$ is defined as the time to reach 90 % of the equilibrium mushy phase thickness $h_e$, and is non-dimensionalized by the diffusive time scale $t^*=1.22\times 10^4\,{\rm s}$. The variations of $t_{e}^{90\,\%}$ with $S_i$ are consistent with those of $h_e$. (e) Thickness $h_e$ shows a clear positive correlation with $t_{e}^{90\,\%}$. A linear fitting gives $h_e/H = 0.0912 \times t_{e}^{90\,\%}/t^*$ (dash-dot line), with correlation coefficient $R^2=0.9846$.

The fact that a more porous ice matrix (high $\phi _e$) can release more salt into the liquid and as a consequence lower the final freezing point of the solution, $T_{0e}$, hence reducing the size of $h_e$, might give a rationale for the globally observed trends. However, this alone is not enough to explain any of the observed fluctuations or non-monotonic dependencies in $h_e$ and $\phi _e$. Moreover, it is important here to notice that the existence of a correlation between $\phi _e$ and $h_e$ is not evident from the present measurements (see figure 3c). If, for instance, mass conservation and geometry were simply ruling the relation between the equilibrium mushy layer thickness and the mushy layer porosity, then one would have $h_e(1-\phi _e) = const.$ (where the constant represents the height of the ice layer in the case of zero salinity and correspondingly no porosity). This implies an increasing trend of $h_e$ with $\phi _e$, while figure 3(c) suggests just the opposite (one can notice a decreasing $h_e(\phi _e)$ for fixed values of $\varTheta _i$). We conclude that the $h_e(\phi _e$) relation is a non-trivial function of the system control parameters $S_i$ and $\varTheta _i$.

Finally, we observe that the trends for $t_{e}^{90\,\%}$ are instead consistent with those of $h_e$ (figure 3d), and a clear linear correlation can be seen between the two quantities (figure 3e). This correspondence can provide a way to estimate the time to equilibrium just by focusing on the magnitude of $h_e$. Therefore, in the following, we will focus on modelling the physical mechanisms responsible for the variations of $h_e$. The question that we wish to address is: can the equilibrium global properties be understood in terms of a model that considers just (i) the role of porosity, and (ii) the role of thermal stratification, on the steady state of the heat transfer process in the system, without taking into account its complex transient dynamics?

5. Multi-layer heat flux model

To answer the above question, we need to understand how $h_e$ and $\phi _e$ are related and how the heat transport across the system is influenced by the convection in the mushy phase and by the combination of stable stratification and unstably stratified turbulence. We will show that a one-dimensional (or integral) multi-layer heat flux model that relies on scaling laws known for natural convection in fluid and porous media can explain $h_e$ with good accuracy once $\phi _e$ is provided as input from the experiment. On the other hand we observe that a quantitative description of $\phi _e$ requires a model for the evolution of the microscopic structure of the porous ice matrix. This is a challenging task that currently goes beyond our capabilities and that we leave for future research. In summary, the multi-layer heat flux model presented in the following allows us to explain the complex relation linking $h_e$ and $\phi _e$, when the latter is considered as an independent variable.

We assume that the position of the solid–fluid interface in the system at the equilibrium $h_e$ is the result of the balance between the mean heat flux across the mushy phase, $\mathcal {F}_m$, and the one in the liquid, $\mathcal {F}_l$ (see figure 4). Both these quantities require a parametrization, which we detail in the following. For conciseness and better readability, we refer the reader to table 1 in Appendix D for the denomination of the many physical quantities and material properties involved in the model.

Figure 4. Synoptic representation of the averaged multi-layer heat flux model. Row 1 provides the adopted parametrization for the global heat flux expression in the mushy phase, $\mathcal {F}_m$, together with the definitions of $Nu_m$ and $Ra_m$. Rows 2–4 report the expression for the global heat flux in the liquid phase, $\mathcal {F}_l$, together with the $Nu_l$ and $Ra_l$ definitions (column 4). Three different combinations of stratifications may exist: (a) unstable, (b) mixed, (c) unstable. The discriminating criterion for their occurrence is the monotonicity of liquid density within the temperature range $[T_{0e}, T_b]$, i.e. the relative magnitude of $T_{0e}$, $T_b$ and the maximum-density temperature $T_{max}$ (column 2), which depends on the salinity $S_e$. Column 3 gives the sketches of these combinations, together with the layer thicknesses and the interface temperatures. The colour gradient indicates the density gradient in the liquid phase, with red representing lighter liquid, and blue representing heavier liquid. The model is supplemented with five input parameters, $T_t$, $T_b$, $S_i$, $\phi _e$ and $H$, and leads to an estimation of $h_e$ (row 5) through the assumption of a statistically steady heat flux balance among layers.

Table 1. Symbols used in this paper.

We express the heat flux across the mushy phase in the statistically steady condition in the form

(5.1)\begin{equation} \mathcal{F}_m = Nu_m \, \varLambda _m(T_m)\,\frac{T_t-T_{0e}}{h_e}. \end{equation}

Equation (5.1) involves the mushy phase Nusselt number $Nu_m$, i.e. the ratio between the total and the conductive heat transfer across the mushy phase. Its dependence as a function of the Rayleigh number $Ra_m$ is known from previous studies on convection in porous media. Here, we use the parametrization by Gupta & Joseph (1973):

(5.2)\begin{equation} Nu_m=\left\{ \begin{array}{@{}ll} 1 , & Ra_m < Ra_{mcr}, \\[2pt] 1.3338+0.0099\,Ra_m, & Ra_m \ge Ra_{mcr}, \end{array} \right . \end{equation}

where the mushy-critical Rayleigh number is $Ra_{mcr}=4{\rm \pi} ^2$. The Rayleigh number for the flow in a porous medium is defined as

(5.3)\begin{equation} Ra_m=\frac{ \varPi \it}{\phi_e}\,\frac{g( \rho(S_t, T_t)- \rho(S_e, T_{0e}))h_e}{\nu(T_m)\,\kappa(T_m)\, \rho(S_e, T_m)}. \end{equation}

This definition involves the evaluation of the permeability $\varPi$, a quantity describing the conductivity of a porous medium to fluid flow (Pal, Joyce & Fleming 2006). The permeability is a function of the structure of the porous ice matrix and is notoriously difficult to measure, in particular for ice, because it is highly sensitive to the environmental conditions and to the history of the medium. Several previous studies, including Kawamura et al. (2006), Petrich, Langhorne & Sun (2006) and Polashenski et al. (2017), revealed that there exists a threshold porosity ($\phi _c$) below which the ice behaves as impermeable. Above such a threshold, the permeability is well approximated by a power-law function of the porosity. We adopt here a similar functional parametrization to relate $\varPi$ to $\phi _e$:

(5.4)\begin{equation} \varPi=\left\{ \begin{array}{@{}ll} 0 , & \phi_e < \phi_{cr}, \\[2pt] C(\phi_e - \phi_{cr})^\alpha, & \phi_e \ge \phi_{cr}, \end{array} \right. \end{equation}

where the three adjustable parameters are chosen as $C=7 \times 10^{-8}$ m$^2$, $\alpha =3$ (as at the leading order of the Carmen–Kozeny relation; Wells et al. 2019) and $\phi _{cr}=0.054$ to best fit $h_e$ measured experimentally. Although these values are ad hoc for our study, the modelled $\varPi$ is of the same order of magnitude as the ones reported in sea ice laboratory and field researches (Kawamura et al. 2006; Petrich et al. 2006; Polashenski et al. 2017).

Similarly, the mean heat flux in the liquid phase can be expressed in the approximate form

(5.5)\begin{equation} \mathcal{F}_l = \left\{ \begin{array}{@{}ll} Nu_l \, \varLambda_l(T_l)\,\dfrac{T_{0e}-T_b}{H-h_e}, & T_{max} \le T_{0e}, \\ \varLambda_l(T_s)\,\dfrac{T_{0e}-T_{max}}{h_{max}} = Nu_l \, \varLambda_l(T_u)\,\dfrac{T_{max}-T_b}{H-h_e-h_{max}}, & T_{max}\in (T_{0e},T_b),\\ \varLambda_l(T_l)\,\dfrac{T_{0e}-T_b}{H-h_e}, & T_{max} \ge T_b. \\ \end{array} \right . \end{equation}

Note that the above expression is composed of three distinct cases according to the value, $T_{max}$, of the temperature corresponding to the maximal density of the liquid solution at the equilibrium with respect to $[T_{0e}, T_b]$. This reflects in three distinct cases: (a) unstable stratification, (b) a stable layer above an unstable one, and (c) a single stable layer; see figure 4. The scale $h_{max}$, appearing in case (b), is the height of the stably stratified layer in the intermediate case. It can be determined by equating the two fluxes across the stable and unstable layers of liquid. The dimensionless heat flux in the liquid is accounted for by the Nusselt number $Nu_l$, dependent on the Rayleigh number $Ra_l$, a relation that is relatively well known. Wang et al. (2021c) proposed the parametrization below by fitting the data from their simulations on the freezing of fresh water in an RB convection system, which leads to theoretical predictions agreeing well with their experiment results. The same parametrization is adopted in the present study:

(5.6)\begin{equation} Nu_l =\left\{ \begin{array}{@{}ll} 1 , & Ra_l < Ra_{lcr}, \\ 0.12+0.88\,Ra_l / Ra_{lcr}, & Ra_l \in [Ra_{lcr}, 1.23\,Ra_{lcr}], \\ 0.27(Ra_l -Ra_{lcr})^{0.27}, & Ra_l > 1.23\,Ra_{lcr}, \end{array} \right . \end{equation}

where the critical Rayleigh number for the onset of convection is $Ra_{lcr}=1708$. We note that the definition of the Rayleigh number in the liquid is needed only for the two cases that involve a density unstable stratification (cases (a) and (b) in figure 4). We use

(5.7)\begin{equation} Ra_{l} = \left\{ \begin{array}{@{}ll} \dfrac{g( \rho(S_e, T_{0e})- \rho(S_e, T_b))( H-h_e )^3 }{\nu(T_l)\,\kappa(T_l)\,\rho(S_e, T_l)}, & T_{max}\le T_{0e},\\ \dfrac{g( \rho(S_e, T_{max})- \rho(S_e, T_b))( H-h_e-h_{max}) ^3 }{\nu(T_u)\,\kappa(T_u)\,\rho(S_e, T_u)}, & T_{max}\in (T_{0e},T_b). \end{array} \right . \end{equation}

Finally, the model tuning requires the numerical values of material properties (see Appendix C for the empirical formulas adopted), which are evaluated at the equilibrium salinity $S_e$. Given the facts that in the present system no salinity source exists, that liquid in the mushy phase contains only a small salt fraction, and that salt is mixed in the liquid either vigorously by convection or slowly by diffusion, we can assume that $S_e$ is uniform in the liquid layer. Hence $S_e$ can be related implicitly to $h_e$ and the inputs of the model by the mass conservation of the salt:

(5.8)\begin{equation} S_e\,\rho(S_e, T_{mean})\,(H-h_e(1-\phi_e))=S_i\,\rho(S_i, T_{b})\,H. \end{equation}

It is worth noticing that (5.8) is suitable only for the present system (when applied to other systems, e.g. open systems, a measurement of $S_e$ must be supplied). The evaluation of $Ra_m$ requires also the value of $\rho (S_t, T_t)$. While $T_t$ is known as it is imposed externally, $S_t$ is not. We assume here that the latter takes the value such that $T_t$ is a local freezing point, $T_0(S_t)=T_t$; in other words, we assume local thermodynamical equilibrium within the mushy phase (Wells et al. 2019). Finally, by solving simultaneously the algebraic equations of heat flux balance $\mathcal {F}_m = \mathcal {F}_l$ and the mass conservation of the salt (5.8), using the measured value $\phi _e$, the model leads to the equilibrium mushy phase thickness $h_e$. We refer to Appendix B for a more detailed explanation of how the model is implemented numerically.

The model results are all illustrated in figure 5. Figure 5(a) shows the mean mushy phase thickness $h_e/H$ as a function of $S_i$ together with its comparison with the experiment results. At $\varTheta _i=1/2$ (blue lines and symbols), where ice thickness decreases for increasing salinity, the agreement with the model is fairly good. The model suggests that in this condition, the heat transfer through the mushy phase remains diffusive (mushy diffusion, MD state). Increasing the salinity, the freezing point is lowered, and this slightly reduces the frozen layer thickness. The role of porosity is marginal as here it stays approximately constant with $S_i$ (see figure 3b). For $\varTheta _i=1/3$ and $3/8$ (red and green lines and symbols), the non-monotonic behaviour is also captured by the model.

Figure 5. Model results for ice/mushy phase thickness compared with the experiments, and identification of heat transport modes. MD/MC stand for diffusive/convective heat transfer in the mushy phase; similarly, LD/LC stand for diffusive/convective heat transport in the liquid phase, while LPC stands for diffusive and convective heat transport due to penetrative convection in the liquid phase, i.e. a combination of diffusive and convective heat transport, as a stable stratification lies above an unstable one. (a) The initial superheat values are respectively $\varTheta _i$ values $1/3$ (red), $3/8$ (green) and $1/2$ (blue), and the initial salinity $S_i$ varies from 0 % to $3.5\,\%$. The calculated $h_e$ (lines with open circles) is close to the experimental results (dots). The spatial variations of mushy phase thicknesses from the experiments (shaded areas) are also included, together with results of reduced models ignoring mushy phase convection ($Nu_m \equiv 1$, dash-dot lines) and ignoring liquid stratification (linear density dependence on temperature $\rho (S_e, T)=\rho (S_e, T_b)\,(1-\beta (S_e, T_b)\,(T-T_b))$, dash-double-dot lines). (b) Phase space diagram based on the two effective control parameters $Ra_m$ and $Ra_l$. The third effective control parameter $S_e$ is not presented as an axis, for readability. Instead, different circle sizes are used to distinguish the different $S_e$. There are in total six heat transport modes of the system. Our experiments fall within five cases. Different circle colours are used to distinguish different $\varTheta _i$. The highest salinity ($S_e \gtrsim 2.7\,\%$) is expected to suppress the LPC regime as it removes the density anomaly of the solution. (c–g) Images from the experiment for typical cases corresponding to the five observed heat transport modes.

It is relevant here to also consider two reduced versions of the model, one ignoring the occurrence of convection in the mushy phase (by assuming $Nu_m \equiv 1$, dash-dot lines in figure 5a), the other neglecting the non-monotonic density stratification in the liquid phase (by assuming linear density dependence on temperature, i.e. $\rho (S_e, T)=\rho (S_e, T_b)\,(1-\beta (S_e, T_b)\,(T-T_b))$, where $\beta (S_e, T_b)$ is the liquid thermal expansion coefficient at the bottom plate, dash-double-dot lines). The comparison of the calculated $h_e$ stemming from the two reduced models allows us to interpret the ‘jump’ or sudden increase in the ice thickness as being associated with the occurrence of convection in the mushy phase (mushy convection, MC state) with increased $\phi _e$ (see figure 3b). Convection strengthens the heat transport, promoting a larger mushy phase growth rate, because $\text {d}h/\text {d}t \propto \mathcal {F}_m - \mathcal {F}_l$. The equilibrium mushy layer thickness, as a result of the competitive balance between $\mathcal {F}_m$ and $\mathcal {F}_l$, also increases because of the enhanced heat transfer efficiency in the mushy phase.

However, the mushy phase convention alone is not sufficient to explain the tendency of $h_e$ (see the dash-double-dot lines in figure 5a). Especially large deviations from experiment results can be seen when $S_i$ is low as they fall into LD (liquid diffusion) and LPC (liquid penetrative convection) regimes (see figure 5b). In the LD state, only diffusive heat transport occurs in the liquid phase, either because the whole liquid phase is stably stratified or because the liquid phase is too thin to support convection. In the LPC state, there is always a stable stratification in contact with the mushy phase above the unstably stratified turbulence, also known as penetrative convection. The presence of this stable layer limits the below convective motion and prevents the mushy phase from a direct penetration of the convection. As an effect, the heat transport across the liquid phase is relatively low, enabling the mushy phase to grow thick. When $S_i$ is increased, the temperature differences decrease in the stably stratified region (from $T_{0e}$ to $T_{max}$) and increase in the unstably stratified turbulence (from $T_{max}$ to $T_b$). As a consequence, the convective motion intensifies, becoming even more turbulent, leading to strong heat transport across the liquid phase and so decreasing $h_e$. Similarly, the temperature difference in the unstable stratification is smaller when $\varTheta _i$ is lower. The unstably stratified turbulence is weaker, and the thermal resistance in the stable stratification above plays a more important role in determining the heat transport of the whole liquid phase. This accounts for the relatively larger errors between the dash-double-dot lines (assuming the entire liquid phase to be unstably stratified) and the experiment results in figure 5(a) when $\varTheta _i= 1/3$ and $3/8$, compared to cases with $\varTheta _i= 1/2$ and the same $S_i$. By further increasing $S_i$, the stable stratification layer disappears, and the heat transport in the liquid phase is purely convective (liquid convection, LC state; see figure 5b).

The remaining discrepancies between the model and the experiments are likely to be related to the modelling of the mushy phase, in particular from neglecting its non-homogeneous internal structure and the uncertainties in the parametrization of its permeability. It is also worth mentioning that the heat transport in the LPC regime can also be clearly described by means of two alternative parameters introduced in Wang et al. (2021a). They are the so-called density inversion parameter (which here reads $(T_{max}-T_{0e})/(T_b-T_{0e})$) and the Rayleigh number relative to the whole fluid layer. It is an alternative but equivalent approach with respect to the one adopted here.

In summary, there are six possible heat transport modes in the system, depending on whether the heat transport is purely diffusive (MD) or convective (MC) in the mushy phase, and whether the heat transport in the liquid phase is purely diffusive (liquid diffusion, LD), diffusive and convective in distinct superposed layers (LPC), or convective (LC). The occurrence of these states is determined by the magnitude of the two control parameters $Ra_m$ and $Ra_l$. The boundaries between purely diffusive (MD and LD) and convective (MC, LPC and LC) heat transport are the critical Rayleigh numbers, respectively $Ra_{mcr}=4{\rm \pi} ^2$ in the mushy phase, and $Ra_{lcr}=1708$ in the liquid. On the other hand, the transition from LPC to LC is due to the disappearance of temperature anomaly in the equation of state of the salty water solution. This occurs when $T_{max}< T_{0e}$, or equivalently, roughly when $S_e \gtrsim 2.7$%. The experimental measurements can be mapped into this posterior three-dimensional phase space, identified by the triplet $(Ra_m, Ra_l, S_e)$; see figure 5(b). It is worth noticing that while $S_i$ varies in $[0,3.5]\,\%$, $S_e$ can vary in a range from 0 % to approximately 7 %. We observe that five out of six possible equilibrium state freezing modes (i.e. the $2 \times 3$ combinations of MD, MC states with LD, LPC, LC ones) can be realized in our experiments and identified by our model. Figures 5(c–g) show photos of the interface profile in these typical cases. We may observe that non-flat interfaces are always associated with convective motions in either the liquid or the mushy phase, or both.

6. Conclusions and outlook

We investigated experimentally the complex physical couplings between the freezing of salty water (ice inception, growth and equilibrium) and convective fluid dynamics transport processes. We revealed that the brine convection in the mushy phase and the combination of stable/unstable stratifications, and so conductive/convective heat transfer, are crucial in the determination of an equilibrium state. The mean mush thickness at the equilibrium can be explained accurately by accounting properly for these two mechanisms via a one-dimensional multi-layer heat flux model that builds on the known global heat transport scaling law properties in a Rayleigh–Bénard system in liquid and porous media.

The occurrence of convection inside the mushy matrix enhances the heat transport and increases the advancement of the freezing front. On the other hand, the stable stratification weakens convection and heat transport in the liquid phase, and sustains a longer diffusive growth of the mushy phase, both leading to a higher mushy thickness. Finally, unstably stratified turbulence promotes convection and heat transfer in the liquid phase, leading to a thinner mushy layer. The resulting heat transport mode at the equilibrium state can be categorized into six regimes, depending on whether convective heat transport exists in the mushy region and whether purely diffusive/convective or both regimes of heat transport exist in the liquid phase.

Many open questions remain to be explored. The multi-layer model needs to be validated further and refined in what concerns the modelling of the permeability and evolution of porosity (which is here taken just as an input from experimental measurements). Furthermore, it would be worthwhile to extend this one-dimensional model to include its temporal behaviour, which seems possible given the fact that its growth rate occurs on time scales that are much longer as compared to the ones of fluid convection; in this sense, a quasi-static approximation, as done for pure ice in Wang et al. (2021d), might work. However, the comprehension of the morphological mushy phase interface, the role of brine channels and related desalination calls for more sophisticated spatial models.

The current study provides a first step in the understanding of this complex system. It also lays a foundation for further studies investigating the interactions between the detailed mushy phase structure, interface morphology, flow dynamics and heat transport performance of the system, as well as exploring the influences of other factors (e.g. the inclination of the system, different geometrical aspect ratios or shapes of the container, or the presence of forced next to natural convection). We expect that our findings can lead to a better understanding and global modelling of industrial, geophysical and climatological processes involving flows in liquid and porous media, and phase change in liquid solutions.