1. Introduction
Melting and freezing in the presence of turbulent flows are fundamental processes in both natural and industrial settings, with wide-ranging applications, including the evolution of polar ice (Stanton et al. Reference Stanton, Shaw, Truffer, Corr, Peters, Riverman, Bindschadler, Holland and Anandakrishnan2013; Dutrieux et al. Reference Dutrieux, Stewart, Jenkins, Nicholls, Corr, Rignot and Steffen2014; Cenedese & Straneo Reference Cenedese and Straneo2023), magma oceans (Labrosse, Hernlund & Coltice Reference Labrosse, Hernlund and Coltice2007; Alboussiere, Deguen & Melzani Reference Alboussière, Deguen and Melzani2010), phase-change materials (Jegadheeswaran & Pohekar Reference Jegadheeswaran and Pohekar2009; Hu et al. Reference Hu, Li, Xu and Fan2022; Choure, Alam & Kumar Reference Choure, Alam and Kumar2023), and the behaviour of icy satellites (Spencer et al. Reference Spencer, Pearl, Segura, Flasar, Mamoutkine, Romani, Buratti, Hendrix, Spilker and Lopes2006; Kang & Flierl Reference Kang and Flierl2020). Accurately estimating melt and freeze rates of ice is crucial for understanding and predicting changes in polar ice, which has direct implications for global climate dynamics (Manabe & Stouffer Reference Manabe and Stouffer1995; Ghil & Lucarini Reference Ghil and Lucarini2020), geoscience (Magorrian & Wells Reference Magorrian and Wells2016; Moon et al. Reference Moon, Sutherland, Carroll, Felikson, Kehrl and Straneo2018) and astrophysics (Butcher et al. Reference Butcher, Balme, Gallagher, Arnold, Conway, Hagermann and Lewis2017; Kang et al. Reference Kang, Mittal, Bire, Campin and Marshall2022). However, the two-way interaction between the moving solid–liquid boundary and the surrounding turbulent flow makes these processes highly complex and difficult to model. Extensive research has explored how various types of turbulent flows interact with solid surfaces during both melting and freezing processes (Davis, Müller & Dietsche Reference Davis, Müller and Dietsche1984; Dietsche & Müller Reference Dietsche and Müller1985; Esfahani et al. Reference Esfahani, Hirata, Berti and Calzavarini2018; Favier, Purseed & Duchemin Reference Favier, Purseed and Duchemin2019; Couston et al. Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2021; Wang et al. Reference Wang, Calzavarini, Sun and Toschi2021b ; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023b ). Studies have focused on interface patterns, equilibrium states, and the intricate dynamics of phase boundaries under freezing and melting conditions, as also summarised in the recent comprehensive review by Du, Calzavarini & Sun (Reference Du, Calzavarini and Sun2024).
In realistic environments, melting and freezing often alternate depending on external conditions, such as daily and seasonal solar cycles that affect sea ice (Perovich & Polashenski Reference Perovich and Polashenski2012), or periodic heating and cooling cycles that influence phase-change materials (Choure et al. Reference Choure, Alam and Kumar2023; Yang et al. Reference Yang, Chong, Liu, Verzicco and Lohse2024) and asymmetric external forces that regulate phase change thermal management (Hu & Fan Reference Hu and Fan2023; Hu et al. Reference Hu, Fan, Gao and Stone2025). This alternation between melting and freezing raises intriguing questions regarding symmetry and reversibility. As a fundamental Stefan problem (Rubinstein Reference Rubinstein1971), one might expect phase transformations to exhibit symmetry – melting as the inverse process of freezing – resulting in a final state that is independent of the initial conditions (Ristroph Reference Ristroph2018). However, ambient flows can break this symmetry, giving rise to complex and even bistable equilibrium states. For example, solar-heated melt ponds in polar regions can exhibit bistability, where the final pond state depends on its initial depth (Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023a ).
Similarly, in the system of Rayleigh–Bénard (RB) convection (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012; Lohse & Shishkina Reference Lohse and Shishkina2023, Reference Lohse and Shishkina2024) – a configuration in which a fluid is heated from below and cooled from above – melting and freezing processes can also yield multiple equilibrium states (Purseed et al. Reference Purseed, Favier, Duchemin and Hester2020). In cases where the melting temperature lies between the temperatures of the top and bottom plates, numerical studies have shown that the interface may settle into one of two distinct states: a diffusive state or a convective state, depending on initial conditions (Purseed et al. Reference Purseed, Favier, Duchemin and Hester2020). Based on this, one might hypothesise that this bistability would disappear in highly turbulent flows, leaving only a single convective equilibrium state. However, earlier experimental studies observed two distinct convective states with different interface patterns during melting and freezing, suggesting a persistent influence of initial conditions even in turbulent regimes (Dietsche & Müller Reference Dietsche and Müller1985). Recent experiments by Wang et al. (Reference Wang, Calzavarini and Sun2021a ) have further revealed that the morphology of ice in an RB set-up depends sensitively on the aspect ratio of the container, underscoring the possibility of multiple stable states under identical external conditions. These findings highlight the complexity of the underlying mechanisms and the need for deeper investigation. Understanding the physical mechanisms behind these observations is critical, as they have implications not only for geophysical and industrial applications, but also for fundamental questions about symmetry-breaking and pattern formation in fluid systems. Further research in this area could improve predictive models of phase-change dynamics in turbulent environments.
In this paper, we present direct numerical simulations (DNS) of melting and freezing in the presence of turbulent convection, aiming to elucidate the asymmetry between these two phase-change processes induced by convective flow. We mainly focus on the equilibrium interface height for melting and freezing, which is an important physical property of the system since it is linked to the heat fluxes by considering the heat flux balance across the solid–liquid interface. In particular, we want to understand the dependence on the control parameters. Our work bridges previous experimental findings with new numerical possibilities and insights, revealing the asymmetry of morphology, equilibrium interface heights, and that corresponding heat flux between melting and freezing indeed persists over a broad range of control parameters. We show that the presence of convective flow disrupts the symmetry, and both qualitative and quantitative analyses suggest that this asymmetry arises from the history-dependent behaviour of the convective flow. The existence of multiple stable states indicates that the exact configuration of the solid–liquid interface is influenced not only by the surrounding environmental conditions, but also by its evolution history. Additionally, we propose a physical model that captures this asymmetry, providing predictions for the equilibrium interface height in both melting and freezing processes. This model advances our understanding of the dynamic interaction between turbulent flows and phase-change interfaces.
The rest of the paper is organised as follows. The set-up and the DNS are introduced in § 2. The numerical results on the asymmetric flow structure and the equilibrium states are discussed in § 3. A theoretical model is provided in § 4. Finally, conclusions and an outlook are presented in § 5.
2. Set-up and control parameters
The simulations are carried out using the Navier–Stokes equations within the Boussinesq approximation and the advection–diffusion equation and heat transfer, and the phase-field equation for the phase transition, commonly used to numerically investigate melting and freezing processes in recent studies (Favier et al. Reference Favier, Purseed and Duchemin2019; Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023c
; Du et al. Reference Du, Calzavarini and Sun2024). For
$Pr\geqslant 1$
, as is relevant here, there are very close similarities between two-dimensional (2-D) and three-dimensional (3-D) RB flows (Van Der Poel et al. Reference van der Poel, Stevens and Lohse2013), therefore we will perform most simulations in two dimensions, also in view of our prior work (Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023b
,
Reference Yang, Howland, Liu, Verzicco and Lohsec
; Xu et al. Reference Xu, Bootsma, Verzicco, Lohse and Huisman2024), which shows quantitative agreement between 2-D and 3-D simulations. We also again conducted some 3-D simulations to compare with analogous 2-D cases, finding qualitative analogies, but also some quantitative differences.
The dimensionless governing equations, constrained by the incompressibility condition
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0$
, are
\begin{align}&\qquad\qquad \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = -\boldsymbol{\nabla }p + \sqrt {\frac {Pr}{Ra}}\,{\nabla} ^2 \boldsymbol{u} + \theta \boldsymbol{e}_z, \end{align}
\begin{align}&\qquad\qquad\,\, \frac {\partial \theta }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\theta =\frac {1}{\sqrt {Ra\,Pr}}\,{\nabla} ^2\theta +St\,\frac {\partial \phi }{\partial t}, \end{align}
\begin{align}& \frac {\partial \phi }{\partial t}=\frac {6}{5\,St\,C\,\sqrt {Ra\,Pr}}\left [{\nabla} ^2\phi -\frac {1}{\epsilon ^2}\phi (1-\phi )(1-2\phi +C\theta )\right ] . \\[12pt] \nonumber \end{align}
Here, the temperature field
$\theta = (T-T_m)/\varDelta$
has been non-dimensionalised using the temperature difference
$\varDelta$
between the bottom plate and the melting temperature
$T_m$
, and the velocity field
$\boldsymbol{u}$
has been non-dimensionalised with the free-fall velocity
$U_{f} = \sqrt {g\beta H \Delta }$
, where
$\beta$
is the thermal expansion coefficient,
$\boldsymbol{e}_z$
is the direction opposite to gravity, and the velocity in the solid phase (
$\phi \gt 0.5$
) is forced to be
$0$
. All lengths have been made non-dimensional with the plate separation
$H$
. The physical control parameters in the equations are the Rayleigh, Prandtl and Stefan numbers, domain aspect ratio, and cooling temperature at the top plate, defined as
\begin{align} Ra=\frac {\beta g \Delta H^3}{\nu \kappa } ,\quad Pr=\frac {\nu }{\kappa } ,\quad St=\frac {\mathcal{L}}{c_p\Delta } ,\quad \varGamma =\frac {W}{H} ,\quad \theta _{c}=\frac {T_{c}}{\varDelta } . \end{align}
Here,
$\nu$
is the kinematic viscosity of the liquid,
$\kappa$
is the thermal diffusivity (which is assumed to be equal in the solid and the liquid),
$g$
is the gravitational acceleration,
$T_c$
the absolute value of the cooling temperature at the top plate,
$c_p$
is the specific heat,
$\mathcal{L}$
is the latent heat, and
$W$
and
$H$
are the horizontal and vertical lengths of the domain.
Figure 1. An illustration of the set-up and parameters with two different initial conditions: (a) melting case with initially mostly solid, which melts from the bottom; (b) freezing case with initially mostly liquid, which freezes from the top.
The response parameters are the Nusselt number
$\textit{Nu}=QH/(K\Delta )$
and the Reynolds number
$Re=UH/\nu$
, which indicate the dimensionless heat transport and flow strength in the system. Here,
$Q$
is the heat flux crossing the system,
$K=\kappa \rho c_p$
is the thermal conductivity in the liquid phase, and
$U=\sqrt {\langle u^2\rangle _V}$
is the time and volume averaged root mean square (r.m.s.) velocity. Note that the thermal diffusivities and conductivities in liquid and solid are the same in our simulations. We also quantify the interface morphology by the spatial-averaged interface height
$\overline {h}$
and the wavenumber of the cellular structure at the interface
$k$
.
The simulation set-up is illustrated in figure 1(a,b). We apply two different initial conditions to explore the melting and freezing dynamics: a predominantly solid phase for the melting case (initially 95 % solid, as shown in figure 1
a), and a predominantly liquid phase for the freezing case (initially 95 % liquid, shown in figure 1
b). The temperature boundary conditions are set to a fixed temperature
$\theta =1$
at the bottom plate,
$\theta =-\theta _{c}$
at the top plate, and
$\theta =0$
for the melting/freezing temperature. When the top temperature is below zero, for the melting case, the solid phase gradually melts upwards from the bottom until equilibrium is reached, leaving a solid layer of a certain thickness above the liquid. For the freezing case, the liquid phase solidifies from the top until it similarly reaches an equilibrium state. For the melting case, we set the initial temperature field in the liquid to linear with random noise perturbation, and different initial settings do not affect the results much due to the small proportion of the liquid phase. For the freezing case, we set the initial temperature field as
\begin{align} \theta &=1-z / h_{\textit{init}}+1 / 2 \sin (2 \pi k_{\textit{init}} / \varGamma x) \sin (\pi z / h_{\textit{init}})^2 \notag\\ &\quad -\theta _c (z-h_{\textit{init}}) / (1-h_{\textit{init}}), \end{align}
where the initial height is set to
$h_{\textit{init}}=0.95$
, and the initial wavenumber is set to
$k_{\textit{init}}=6$
, a typical configuration on a Fourier basis that allows the flow and morphological structures to be sustained from the initial state, commonly used to control the initial number of rolls (Wang et al. Reference Wang, Verzicco, Lohse and Shishkina2020). We also further examined the effect of varying
$k_{\textit{init}}$
, and found that it can lead to different final equilibrium states.
We applied periodic boundary conditions in the horizontal direction. For velocity fields, we apply no-slip at the plates. The parameter range that we will explore is for Rayleigh numbers
$10^7\leqslant Ra\leqslant 3\times 10^8$
and for temperatures at the top plate
$-40\leqslant \theta _{t}\leqslant 0$
. By default,
${\textit{Pr}}$
is set to
$1$
– unless specified when we investigate the effect of
${\textit{Pr}}$
and vary in
$1\leqslant Pr\leqslant 10$
, representative of properties of the different fluid –
$\varGamma$
is fixed to
$8$
, and
$St$
is fixed to
$1$
. We also checked the effect of
$St$
for the asymmetry behaviour in Appendix A.
The DNS were performed using a finite difference code (Verzicco & Orlandi Reference Verzicco and Orlandi1996) with the criteria for the grid resolution following Ostilla-Mónico et al. (Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015). The code has been tested extensively, and benchmarked against other codes (Kooij et al. Reference Kooij2018) and applied in various fluid systems (Yang, Verzicco & Lohse Reference Yang, Verzicco and Lohse2016; Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018; Chong et al. Reference Chong, Ng, Hori, Yang, Verzicco and Lohse2021; Yang et al. Reference Yang, Ng, Chong, Verzicco and Lohse2022). We implement the phase-field method (Hester et al. Reference Hester, Couston, Favier, Burns and Vasil2020) to simulate the melting/freezing solid (Liu et al. Reference Liu, Ng, Chong, Lohse and Verzicco2021). In this method, the phase-field variable
$\phi$
is integrated in time and space, and smoothly transitions from value
$1$
in the solid to value
$0$
in the liquid. The applied phase-field model was initially derived based on entropy conservation, which guarantees the thermodynamic consistency and also satisfies the Gibbs–Thompson relation.
Figure 2. Instantaneous snapshots of the temperature field and the solid–liquid interface contour from 2-D simulations of (a) the melting case, and (b) the freezing case at
$Ra=10^7,\ \theta _c=8$
, and the temperature field and solid–liquid interface contour from 3-D simulations of (c) the melting case and (d) the freezing case at
$Ra=10^7,\ \theta _c=8$
. (e) The time evolution of the mean height
$\overline {h}$
for the cases in (a–d).
3. Flow structures and equilibrium states for melting and freezing
3.1. Asymmetric flow structures and equilibrium height
In figures 2(a,b), we present the flow structures and phase boundaries for both the melting and freezing cases under identical control parameters:
$Ra=10^7,\ \theta _c=8$
, as obtained from 2-D simulations. A clear distinction emerges in the number of thermal plumes and the corresponding cellular patterns at the interface. This difference in flow and interface structure between melting and freezing is not limited to 2-D simulations; similar variations are also observed in 3-D simulations with
$Ra=10^7,\ \theta _c=8$
, as shown figure 2(c,d).
These differences extend beyond just flow patterns and interface wavenumbers. The final equilibrium heights
$\overline {h}^m_{\infty }$
and
$\overline {h}^f_{\infty }$
of the phase boundary also differ between melting and freezing processes. Figure 2(e) illustrates the temporal evolution of the spatially averaged interface height
$\overline {h}(t)$
for both melting (
$\overline {h}^m_{\infty }$
) and freezing (
$\overline {h}^f_{\infty }$
) in the cases shown in figure 2(a–d). In an idealised Stefan problem with no ambient flow, one would expect both processes to reach an identical final state, controlled by the balance of heat fluxes between the solid and liquid phases. However, with convective flow, our results reveal two distinct equilibrium heights,
$\overline {h}^m_{\infty }$
and
$\overline {h}^f_{\infty }$
, for melting and freezing, respectively. Interestingly, in the melting case,
$\overline {h}^m(t)$
continues to increase over time, eventually surpassing the steadily decreasing
$\overline {h}^f(t)$
observed in the freezing case, so that
$\overline {h}^m_{\infty }\gt \overline {h}^f_{\infty }$
.
This difference in the equilibrium heights
$\overline {h}^m_{\infty }$
and
$\overline {h}^f_{\infty }$
highlights an asymmetry induced by the convective flow, which is a marked difference from traditional Stefan problem expectations. The presence of buoyancy-driven flow not only modifies the flow-field structure and interface geometry, but also influences the overall phase-change dynamics, leading to different thermal equilibria. This finding underscores the critical role of convective effects in breaking the symmetry between melting and freezing, suggesting that the history and intensity of convective flows can have lasting impacts on the final interface configuration, which is relevant to a variety of natural and industrial applications.
3.2. Asymmetry dependence on
$Ra$
and initial condition
Figure 3. (a) The mean height
$\overline {h}_{\infty }$
as a function of
$\theta _c\, Ra^{-1/4}$
for melting and freezing cases with different
$Ra$
. The dashed line represents the location when the curves for freezing cases transition to increase. (b) The wavenumber
$k/k_0$
as a function of
$\theta _c\, Ra^{-1/4}$
for melting and freezing cases with different
$Ra$
, where
$k_0$
is the wavenumber for
$\theta _c=0$
. (c) The r.m.s. of the height,
$h_{rms,\infty }$
, as a function of
$\theta _c\, Ra^{-1/4}$
for melting and freezing cases with different
$Ra$
. The dashed vertical lines in (a–c) are at the same value,
$\theta _c\,Ra^{-1/4}=0.2$
.
The observed asymmetry between melting and freezing is evident across a broad range of Rayleigh numbers
$Ra$
and cooling temperatures
$\theta _c$
. In figure 3(a), we plot the equilibrium mean interface heights
$\overline {h}^m_{\infty }$
and
$\overline {h}^f_{\infty }$
as functions of the rescaled parameter
$\theta _c\, Ra^{-1/4}$
across various values of
$Ra$
and
$\theta _c$
. This rescaling is based on the heat flux balance at the interface, as well as the classical heat transfer scaling
$\textit{Nu}\sim Ra^{1/4}$
from Grossmann–Lohse theory (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001), which compensates for the effect of
$Ra$
, and facilitates a clearer comparison between cases. Consistently, the melting cases exhibit higher values of
$\overline {h}_{\infty }$
than the freezing cases under the same conditions of
$Ra$
and
$\theta _c$
, i.e.
$\overline {h}^m_{\infty }\gt \overline {h}^f_{\infty }$
. Notably, the difference in
$\overline {h}_{\infty }$
between the melting and freezing cases increases with larger values of
$\theta _c$
, and decreases with decreasing
$\overline {h}_{\infty }$
. Beyond a certain threshold of
$\theta _c$
, however, this difference begins to diminish. This is because when the liquid layer becomes sufficiently shallow, and the flow approaches a laminar state, the differences between the melting and freezing cases eventually vanish. At higher
$Ra$
, the effective
$Ra$
remains larger for the same
$\overline {h}_{\infty }$
, meaning that the differences persist until
$\overline {h}_{\infty }$
becomes even smaller, or
$\theta _c$
becomes even larger. We will analyse this in further detail in the following subsections.
Figure 3(b) presents the wavenumber
$k$
of the cellular structure at the interface as a function of
$\theta _c\, Ra^{-1/4}$
. Here again, the melting cases show consistently higher wavenumbers than the freezing cases for identical values of
$Ra$
and
$\theta _c$
, with
$k$
increasing as
$\theta _c$
rises. This correlation between
$\overline {h}^{m,f}_{\infty }$
and
$k$
reflects a broader trend: systems with a larger number of convection rolls exhibit enhanced heat transport properties, as seen in prior studies on RB convection (Wang et al. Reference Wang, Verzicco, Lohse and Shishkina2020) and Taylor–Couette flow (Huisman et al. Reference Huisman, van der Veen, Sun and Lohse2014). Based on the heat flux balance between the liquid and solid phases, improved heat transport within the liquid phase translates to a higher equilibrium height
$\overline {h}^{m}_{\infty }$
during the melting process.
Note that the statistics for the freezing cases in figure 3(b) appear to be more chaotic at higher cooling temperatures. This is because, as the cooling temperature increases, splitting events occur intermittently, allowing
$k$
to remain approximately constant within a certain range of
$\theta _c$
until the next splitting event occurs. Furthermore, at higher cooling temperatures, the fluid layer becomes shallower, and the flow weakens, resulting in longer time scales for new splitting events to develop.
Figure 3(c) shows the morphology roughness, represented by the r.m.s. of the height,
$h_{rms,\infty }$
, as a function of
$\theta _c\, Ra^{-1/4}$
. The freezing cases show consistently rougher interfaces than the melting cases, given the same control parameters. Meanwhile, for the freezing cases,
$h_{rms,\infty }$
shows a non-monotonic trend. As
$\theta _c$
increases,
$h_{rms,\infty }$
first increases and then decreases, with the same transition location as the increase of
$k$
. Note that there is a difference between the 2-D and 3-D results. The difference between freezing and melting is relatively smaller in three dimensions than in two, likely due to flow patterns in the third direction destabilising the 2-D circulation, leading to merging and splitting events. This effect reduces the observed wavenumber difference between melting and freezing cases. However, we believe that the difference in equilibrium heights between freezing and melting remains robust, as it has also been observed in experiments. For instance, Dietsche & Müller (Reference Dietsche and Müller1985) reported
$12$
cells in the freezing case, and
$15$
cells in the melting case, under identical control parameters. Similarly, Wang et al. (Reference Wang, Calzavarini and Sun2021a
) demonstrated that multiple stable ice profiles can emerge, corresponding to different convective states under the same external conditions. Furthermore, we will present a model based on heat flux balance to quantitatively unify all results in § 4, including both 2-D and 3-D results.
To explore the role of the initial flow structure in the freezing cases, we performed additional tests by initialising the freezing cases with varying numbers of rolls, each corresponding to a different initial wavenumber
$k_{\textit{init}}$
. These cases demonstrate that the initial morphology can indeed be sustained, with the number of rolls persisting through the evolution, as illustrated in figure 4(a–d). In figure 4(e), we plot
$\overline {h}^{f}_{\infty }$
as function of
$\theta _c\, Ra^{-1/4}$
for different initial roll numbers
$k_{\textit{init}}$
, revealing an intriguing trend: the threshold of
$\theta _c$
for the transition in the freezing curve and the associated
$\overline {h}^{f}_{\infty }$
vary depending on the initial wavenumber, as depicted in the inset in figure 4(e). This resembles the case of standard RB turbulence, for which Wang et al. (Reference Wang, Verzicco, Lohse and Shishkina2020) showed that for a given aspect ratio and
$Ra$
,
${\textit{Pr}}$
, only a certain number of rolls form a stable equilibrium.
Figure 4. (a–d) Instantaneous snapshots of temperature fields with the solid–liquid interface contour for different
$k_{\textit{init}}$
at
$Ra=10^7$
at the equilibrium stage. (e) The mean height
$\overline {h}^{f,m}_{\infty }$
as a function of
$\theta _c\, Ra^{-1/4}$
for freezing cases with different
$k_{\textit{init}}$
, and the melting case, at
$Ra=10^7$
. The dashed lines show the locations where the curves for freezing cases transition and increase. We also include the melting case
$\overline {h}^{m}_{\infty }$
. The inset represents
$\overline {h}^{f}_{\infty }$
as a function of
$k_{\textit{init}}$
. The black points represent the transition points, and the solid line represents the prediction from (3.2).
We further test different initial conditions for the freezing case by setting the flow with various numbers of rolls initially. Given a certain number of rolls
$k_{\textit{init}}$
initially, the morphology can be adjusted and sustained to follow the number of rolls, as shown in figure 4(a–d). Figure 4(e) shows
$\overline {h}^{f,m}_{\infty }$
as a function of
$\theta _c\,Ra^{-1/4}$
with different initial number of rolls. Interestingly, these cases show the transition at different
$\theta _c$
and corresponding
$\overline {h}^{f}_{\infty }$
, which is also shown in the inset.
3.3. Equilibrium states predicted by the critical aspect ratio
$\varGamma _c$
The stability of initial conditions can be quantified by the aspect ratio of the convection cells, denoted
$\varGamma _c={\varGamma /(k\overline {h}_\infty )}$
. For the transient melting processes in RB convection, Yang et al. (Reference Yang, Howland, Liu, Verzicco and Lohse2023c
) showed that
$\varGamma _c$
is bounded within a specific range, constrained by the elliptical instability criteria for convection cells (Wang et al. Reference Wang, Verzicco, Lohse and Shishkina2020; Zwirner, Tilgner & Shishkina Reference Zwirner, Tilgner and Shishkina2020). In brief, the strain–vorticity balance in 2-D flows is determined by the (elliptical) shape of the rolls. The aspect ratio of the rolls is directly related to the strain and the vorticity. To form a roll, the strain must be bounded by the vorticity, which gives the limits of
$\varGamma _c$
. This range provides theoretical lower and upper bounds for the aspect ratio that allow stable convection patterns within the system. Yang et al. (Reference Yang, Howland, Liu, Verzicco and Lohse2023c
) further showed that the merging events in melting cases occur when the aspect ratio decreases until it reaches the lower bounds.
Figure 5. The aspect ratio
$\varGamma _c$
as a function of the effective Rayleigh number
$Ra_{\textit{eff}}$
for all simulation cases.
We examine these bounds in the equilibrium state by analysing the equilibrium structure of the solid–liquid interface for various control parameters. Figure 5 illustrates
$\varGamma _c$
as a function of
$Ra_{\textit{eff}}$
for both melting and freezing cases, where
$Ra_{\textit{eff}}=\alpha g\Delta \overline {h}_{\infty }^3/(\nu \kappa )$
is the effective Rayleigh number, which serves as an essential measure of the effective buoyancy-driven convection strength in the system. The plot reveals two distinct boundaries, with theoretical bounds derived based on the work of Shishkina (Reference Shishkina2021), which gives
\begin{equation} \varGamma _{c1}\lesssim \varGamma _c\lesssim \varGamma _{c2}, \end{equation}
where
$\varGamma _{c1}=0.7$
and
$\varGamma _{c2}=1.6$
are independent of
$Ra_{\textit{eff}}$
. For the same control parameters, the aspect ratios in the melting cases tend to be lower and closer to the lower bound compared to those in the freezing cases. This finding aligns with our observations of distinct equilibrium heights and wavenumber characteristics in the two processes in figure 3. Melting cases favour compact convection cells with lower aspect ratios due to the growth and merging from small plumes, while freezing cases tend to maintain the elongated cells with higher aspect ratios from the initial large plumes, indicative of differing heat transport efficiencies.
This analysis also provides insight into the possible range of initial wavenumbers, as the system’s behaviour is constrained within these theoretical bounds of
$\varGamma _c$
. For initial states falling within the lower and upper bounds, stable roll structures are more likely to form, while initial conditions outside this range lead to unstable or transitional flow patterns.
3.4. Asymmetry dependence on
${\textit{Pr}}$
Another important control parameter in this study is the Prandtl number
${\textit{Pr}}$
, which we have fixed as
$1$
to maintain consistency. However, in real-world applications,
${\textit{Pr}}$
depends on fluid properties that can vary significantly, such as reaching value
$Pr=10$
for cold water near
$0\,^\circ\textrm{C}$
. Figure 6 shows the effect of varying
${\textit{Pr}}$
from
$1$
to
$10$
on the observed asymmetry in melting and freezing behaviours. Notably, as
${\textit{Pr}}$
increases, the asymmetry between melting and freezing becomes less pronounced, and ultimately vanishes at
$Pr=10$
, as depicted in figure 6(d), which means that for ice melting in cold water, this asymmetry should be negligible.
Figure 6. The effect of Prandtl number. The 2-D simulations of melting and freezing cases at
$Ra=10^7$
and (a)
$Pr=1$
, (b)
$Pr=4$
, (c)
$Pr=10$
. (d) The mean height
$\overline {h}_{\infty }$
as a function of
${\textit{Pr}}$
for different
$\theta _c$
, for both melting and freezing.
The changes in flow structure at different
${\textit{Pr}}$
values shed light on this phenomenon. At
$Pr=1$
, as shown in figure 6(a), the flow is more coherent, with a circulation-dominated structure that supports large-scale convective rolls. This organised flow pattern promotes distinct cellular structures at the solid–liquid interface, enhancing the asymmetry between melting and freezing. However, as
${\textit{Pr}}$
increases, the flow becomes more and more chaotic (see figure 6
b,c), with increasing number of plumes, and disrupting the large-scale circulation. Smaller plume-dominated structures emerge, creating a more randomised and chaotic flow. In this scenario, the melting case, which initially has more cells, tends to experience merging events, reducing the total number of cells. Conversely, in the freezing case, where fewer cells are present, splitting events may occur, increasing the number of cells. As a result, the contrast between melting and freezing cases diminishes as the flow self-reorganises dynamically, leading to less asymmetry between the two cases at higher Prandtl numbers.
This shift from coherent to chaotic flow with higher
${\textit{Pr}}$
underscores the influence of the fluid properties on phase-change processes. A high
${\textit{Pr}}$
inhibits the development of stable cellular patterns, resulting in a loss of the organised structure that initially drives the asymmetry. As a result, for applications where
${\textit{Pr}}$
is naturally high, such as in cold-water dynamics, the distinction between melting and freezing dynamics may be less evident. However, it requires a more comprehensive exploration of the parameter space for different
${\textit{Pr}}$
to fully understand the effect of
${\textit{Pr}}$
, including the influence of
${\textit{Pr}}$
at different
$\theta _c$
. In the remainder of this study, we maintain
$Pr=1$
to focus on the regime where asymmetry is most pronounced, and examine the underlying mechanisms in more detail.
3.5. Convection cell splitting events in freezing cases
Another intriguing aspect of this study is the complex evolution observed in the freezing case. While the melting dynamics has been well-documented in the literature, with studies showing that the solid continuously melts as underlying plumes grow, merge and restructure (Favier et al. Reference Favier, Purseed and Duchemin2019; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023c
), the freezing process exhibits distinctive, non-monotonic behaviour, particularly when
$\theta _c$
is sufficiently high. As shown in figure 7(b), the evolution of the mean interface height
$\overline {h}^f(t)$
during freezing displays a non-monotonic trend, indicating a series of dynamic structural transitions that are not present in the melting case. Similar behaviour for
$h_{rms,\infty }$
can be seen in figure 7(c,d).
To gain further insight into this feature, we examine the temporal evolution of the interface height for both melting and freezing cases, presented in figure 7(e, f). At
$\theta _c=7$
, the melting case shows the expected merging of convection cells over time until a stable equilibrium state is achieved. In contrast, the freezing case maintains a stable cellular pattern from start to finish, with no significant changes in cell number or arrangement. However, at a larger top cooling temperature,
$\theta _c=14$
, for the melting case, convection rolls continue to exhibit merging behaviour. Interestingly, for the freezing case, convection rolls undergo a notable ‘splitting’ event, where existing cellular structures divide, forming additional cells. This splitting corresponds to a point in time when
$\overline {h}^f(t)$
begins to increase, indicating a secondary phase in the freezing process.
Figure 7. The temporal evolution of the mean height
$\overline {h}$
for melting and freezing cases at the same control parameters: (a)
$Ra=10^7,\ \theta _c=7$
, and (b)
$Ra=10^7,\ \theta _c=14$
. In both cases,
$Pr=1$
. The dashed line in (b) represents the time when the curve of the freezing case starts to increase. (c,d) The temporal evolution of
$h_{rms}$
for melting and freezing cases at the same control parameters. The temporal evolution of the interface height
$h(x)$
for melting (left) and freezing (right) cases at (e)
$Ra=10^7,\ \theta _c=7$
, and ( f)
$Ra=10^7,\ \theta _c=14$
. The horizontal dashed line in ( f) represents the same time as in (b) and (d) when the curve of the freezing case starts to increase.
The nature of these splitting events becomes clearer when we examine the flow and temperature fields, as illustrated in figure 8. Figure 8(a) depicts the convection roll merging process over time during melting, while figure 8(b) focuses on the convection roll splitting dynamics during freezing. Consistent with prior studies (Favier et al. Reference Favier, Purseed and Duchemin2019; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2023c ), the cellular structures at the interface are supported by the convective plumes and circulatory flow underneath. In the melting case, the stronger heat flux from rising hot plumes compared to sinking cold plumes results in a more pronounced melt rate beneath each convection cell. As melting progresses, plumes can grow and destabilise, leading to the merging of adjacent cells.
In contrast, during the freezing process, a different mechanism emerges when
$\theta _c$
is high enough. As the solid phase advances towards the bottom plate, it creates a thin, weak-flow region close to the solid–liquid boundary. This region provides a low-shear environment that helps the formation of new plumes. As these new plumes grow, they efficiently enhance heat transport, increasing the local heat flux and leading to further melting, forming a new cellular structure. As the cell continues to melt and grow, it eventually merges into a larger roll, effectively ‘splitting’ the existing cells at the interface. Figure 8(c) illustrates this convection roll splitting phenomenon, which occurs only when
$\theta _c$
reaches a threshold that brings the freezing front into close proximity with the bottom plate. This splitting behaviour in the freezing case underscores a key difference in the underlying mechanisms governing melting and freezing in convective environments. While melting processes predominantly result in the merging of convection cells, freezing introduces conditions that facilitate cell division and increased structural complexity. This difference in behaviour at higher
$\theta _c$
values is a direct outcome of the two-way coupling between the phase boundary and the convective flow.
Figure 8. Zoom-in (height from
$0$
to
$0.5$
) snapshots of the temperature field, the velocity vectors, and the solid–liquid interface contour for (a) melting and (b) freezing cases at
$Ra=10^7,\ \theta _c=14$
. From top to bottom represents the time evolution of merging cells for the melting case, and splitting cells for the freezing case. (c) Illustrations of the splitting cell event for the freezing case.
We further provide a theoretical model to predict the occurrence of convection roll splitting events. The non-monotonic behaviour of the
$\overline {h}^f(t)$
in freezing cases is due to the splitting of convection rolls, which enhances the heat flux, leading to an intermediate increase in
$\overline {h}^f(t)$
. The splitting is caused by the instability of convection rolls. As the solid keeps freezing,
$\overline {h}^f(t)$
decreases, and the instantaneous aspect ratio increases until it reaches the upper limit and cannot be sustained any longer, then the splitting event occurs. Therefore, the splitting of convection rolls is directly related to the stability of these rolls, and can be predicted from the critical aspect ratio
$\varGamma _c$
in (3.1). This critical upper bound explains the occurrence of splitting events in the freezing cases, and gives a prediction for the critical height
$\overline {h}_c$
when the splitting event occurs as a function of the initial wavenumber
$k_{\textit{init}}$
, as
\begin{equation} \overline {h}_c=\frac {\varGamma }{\varGamma _{c2} k_{\textit{init}}}, \end{equation}
which is plotted in the inset of figure 4(b), and shows good agreement with the transition points (corresponding to when the splitting event occurs) in numerical simulations.
4. Unifying model for the equilibrium state
$\overline {\boldsymbol{h}}_{\boldsymbol{\infty}}$
as a function of
$\boldsymbol{\theta} _{\boldsymbol{t}}$
,
$\boldsymbol{Ra}$
and
$\boldsymbol{\varGamma}_ {\boldsymbol{c}}$
Building on the introduced effective Rayleigh number
$Ra_{\textit{eff}}$
and the cell aspect ratio
$\varGamma _c$
, we further propose a unifying model to predict the mean height
$\overline {h}_{\infty }$
under varying
$\theta _{t}$
,
$Ra$
and
$\varGamma _c$
. The equilibrium height can be derived from the heat flux balance across the solid–liquid interface, expressed as
\begin{equation} \textit{Nu}\, \frac {1}{\overline {h}_{\infty }}=\frac {\theta _c}{1-\overline {h}_{\infty }}, \end{equation}
where the term on the left-hand side represents the heat flux in the liquid phase, and the term on the right-hand side corresponds to the heat flux in the solid phase. Here,
$\textit{Nu}$
remains unknown and needs to be determined to solve for
$\overline {h}_{\infty }$
.
In classical RB convection,
$\textit{Nu}$
can be well predicted using the Grossmann–Lohse theory, which predicts
$\textit{Nu}\sim Ra^{1/4}$
for the range of
$Ra$
and
${\textit{Pr}}$
used in our study. Furthermore, to account for the effect of
$\varGamma _c$
, we modify this scaling relationship by assuming that
$\textit{Nu}$
has a linear dependence on
$\varGamma _c$
. We acknowledge that this linear dependence does not hold universally, and refer to Shishkina (Reference Shishkina2021) for a more detailed analysis of the
$\textit{Nu}$
dependence on
$\varGamma _c$
. However, in our parameter space, plotting the rescaled quantity
$\textit{Nu}\,Ra^{-1/4}$
against
$\varGamma _c$
reveals a linear trend, and the same for
$Re\,Ra^{-1/2}$
against
$\varGamma _c$
, as shown in the insets of figure 9(a,b), allowing us to approximate
$\textit{Nu}$
as
\begin{align}& {\textit{Nu}} =a_0\,Ra^{1/4}_{\textit{eff}}\,Pr^{-1/12}\,f_{\textit{Nu}}(\varGamma _c), \end{align}
\begin{align}&\,\,\, Re =b_0\,Ra^{1/2}_{\textit{eff}}\,Pr^{-1}\,f_{Re}(\varGamma _c), \end{align}
where
$f_{\textit{Nu}}(\varGamma _c)=1-a_1(\varGamma _c-\varGamma _{c0})$
and
$f_{Re}(\varGamma _c)=1-b_1(\varGamma _c-\varGamma _{c0})$
are linear functions of
$\varGamma _c$
, with
$\varGamma _{c0}$
fixed as
$1$
, corresponding to the typical aspect ratio for the Grossmann–Lohse theory. The constants
$a_0$
and
$b_0$
are empirically determined coefficients, and
$a_1$
and
$b_1$
represent the sensitivity of
$\textit{Nu}$
and
$Re$
to variations in
$\varGamma _c$
, which are obtained by the fitting of the simulation data.
Figure 9. (a) Plot of
$\textit{Nu}\,f^{-1}_{\textit{Nu}}$
as a function of
$Ra_{\textit{eff}}$
from various numerical simulations; see legend in (c). The solid line shows the scaling relation of (4.2), and the inset shows
$\textit{Nu}\,Ra^{-1/4}$
against
$\varGamma _c$
, with the solid line being the linear fitting. (b) Plot of
$Re\,f^{-1}_{Re}$
as a function of
$Ra_{\textit{eff}}$
. The solid line is the scaling relation (4.3), and the inset shows
$Re\,Ra^{-1/2}$
against
$\varGamma _c$
, with the solid line being the linear fitting. (c) The equilibrium interface height
$\overline {h}_{\infty }$
as a function of
$\theta _c\,\textit{Nu}_0^{-1}\,f_{\textit{Nu}}^{-1}(\varGamma _c)$
based on (4.4).
To validate this model for the
$\varGamma _c$
dependence, we plot the compensated quantities
$\textit{Nu}\,f^{-1}_{\textit{Nu}}$
and
$Re\,f^{-1}_{Re}$
as functions of
$Ra_{\textit{eff}}$
in figure 9(a,b). The results show that across a range of cases, the data collapse onto a single scaling curve, as predicted by (4.2) and (4.3), which supports the robustness of the linear relation between
$\textit{Nu}$
,
$Re$
and
$\varGamma _c$
for the parameter space in our study.
By combining (4.1) and (4.2), we derive an expression for
$\overline {h}_{\infty }$
:
\begin{equation} \overline {h}_{\infty }^{-1/4}(1-\overline {h}_{\infty })=\theta _c\,\textit{Nu}_0^{-1}\,f_{\textit{Nu}}^{-1}(\varGamma _c), \end{equation}
where
$\textit{Nu}_0=a_0\,Ra_{\textit{eff}}^{1/4}$
. In figure 9(c), we plot
$\overline {h}_{\infty }$
as a function of the rescaled cooling temperature
$\theta _c\,\textit{Nu}_0^{-1}\,f_{\textit{Nu}}^{-1}(\varGamma _c)$
. The black line represents the theoretical prediction from (4.1), showing excellent agreement with the simulation data.
5. Conclusion and outlook
In conclusion, this study provides a detailed investigation into the asymmetric behaviours of melting and freezing processes in RB convection, revealing distinct differences in the equilibrium flow structures, solid–liquid interface morphology, and the equilibrium mean interface heights
$\overline {h}^m_{\infty }$
and
$\overline {h}^f_{\infty }$
. Our results demonstrate that this asymmetry is robust across a wide range of Rayleigh numbers
$Ra$
and top cooling temperatures
$\theta _c$
, indicating that the convective dynamics plays a critical role in breaking the symmetry between these phase-change processes. A close examination of the morphological evolution in the freezing cases reveals the occurrence of a convection roll splitting event in the cellular structures at the interface, which accounts for the observed non-monotonic trend in the height evolution. This convection roll splitting, absent in the melting cases, provides insights into the role of convective instabilities and plume dynamics in shaping the solid–liquid boundary during freezing. Such events lead to the formation of additional convective cells, and explain the sustained difference in interface height between the two processes.
To quantify these morphological differences, we have introduced an effective aspect ratio
$\varGamma _c$
of the cellular structures at the interface, which is shown to be constrained between theoretical upper and lower bounds derived from the convection roll stability criteria of Wang et al. (Reference Wang, Verzicco, Lohse and Shishkina2020) and Shishkina (Reference Shishkina2021). The observed aspect ratio bounds align with theoretical predictions, reinforcing the idea that the convective cell structure in the liquid phase strongly influences the overall interface morphology. Finally, using the Grossmann–Lohse theory and the heat flux balance, we develop a unifying model that accurately predicts the equilibrium interface height
$\overline {h}$
across various initial conditions and control parameters. The resulting predictions are consistent with our numerical results, successfully capturing the trends in
$\overline {h}_{\infty }$
across melting and freezing cases,
$\overline {h}^m_{\infty }$
and
$\overline {h}^f_{\infty }$
.
These findings advance our understanding of the complex interactions between phase-change dynamics and convective flow, highlighting the importance of historical flow states and morphological evolution in determining the final equilibrium conditions. This work not only extends existing models of convective phase change by accounting for asymmetry in melting and freezing, but also provides a framework for exploring similar behaviours in more complex, 3-D systems and applications involving phase-change materials and thermal management technologies.
Future studies could explore how varying boundary conditions (Yang et al. Reference Yang, Chong, Liu, Verzicco and Lohse2024), initial conditions (Purseed et al. Reference Purseed, Favier, Duchemin and Hester2020; Wang et al. Reference Wang, Calzavarini and Sun2021a
), external forces (Couston et al. Reference Couston, Hester, Favier, Taylor, Holland and Jenkins2021) and saline concentration (Du et al. Reference Du, Wang, Jiang, Calzavarini and Sun2023; Xue, Zhang & Ni Reference Xue, Zhang and Ni2024), as well as freshwater temperature with density anomaly effect (Wang et al. Reference Wang, Calzavarini, Sun and Toschi2021b
), further influence the symmetry-breaking mechanisms in convective melting and freezing, as well as the merging and splitting events. The persistence of asymmetry at higher
$Ra$
also warrants further investigation, as the increasingly chaotic plume and circulation patterns may diminish the differences between the melting and freezing cases, similar to the observed trends at high
${\textit{Pr}}$
.
Funding
This work was financially supported by the European Union (ERC, MultiMelt, 101094492). It was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. We also acknowledge the EuroHPC Joint Undertaking for awarding the project EHPC-REG-2023R03-178 access to the EuroHPC supercomputer Discoverer.
Declaration of interests
The authors report no conflict of interest.
Appendix A. The effect of
$\boldsymbol{St}$
on the asymmetric melting and freezing
Considering that the Stefan number
$St$
is different in realistic cases as we have modelled here (
$St=1$
for computational efficiency), we conducted additional simulations with
$St=10$
(e.g. ice melting in
$8\,^\circ \textrm {C}$
water) and
$St=4$
(e.g. ice melting in
$ 20\,^\circ \textrm {C}$
water) at
$Ra=10^7$
,
$Pr=1$
and
$\theta _c=8$
to examine the influence of
$St$
on the melting and freezing processes. Figure 10 presents the temporal evolution of the spatially averaged interface height
$\overline {h}(t)$
for these cases. As anticipated, the time required to reach equilibrium increases with a larger
$St$
, reflecting the larger energy needed for phase change under higher
$St$
conditions. This energy requirement leads to a slower rise of
$\overline {h}_\infty$
during melting, which, in turn, allows more time for plume merging and reorganisation, ultimately resulting in a lower equilibrium wavenumber
$k$
and a slightly reduced final
$\overline {h}_\infty$
.
Figure 10. Temporal evolution of
$\overline {h}(t)$
at
$Ra=10^7$
,
$Pr=1$
and
$\theta _c=8$
for different
$St$
.
Despite these effects on the dynamics, the asymmetry between melting and freezing processes persists across varying
$St$
. Both our observations and the theoretical framework that we developed remain consistent with changes in
$St$
. Thus although
$St$
affects the temporal evolution and the time scale of phase change, it does not significantly alter the underlying asymmetry in interface morphology or the flow structure between melting and freezing cases. Given this result, we have chosen to set
$St=1$
in our simulations, which allows for faster computational convergence without compromising the generality of our findings.
Appendix B. Resolution test for the evolution of mean height
$\overline {\boldsymbol{h}}$
To ensure the adequacy of our chosen resolution, we perform simulations at different resolutions to assess convergence, as shown in figure 11. We select
$Ra=10^7$
and
$Pr=1$
with two different
$\theta _c$
, as
$\theta _c=7$
(without splitting event) and
$\theta _c=14$
(with splitting event). The results demonstrate the convergence of the temporal curve as the resolution increases. Based on this, we finalise our choice of vertical resolution as
$192$
, constant for different
${\textit{Pr}}$
, and doubling it for
$Ra=10^8$
.
Figure 11. Temporal evolution of
$\overline {h}(t)$
at
$Ra=10^7$
,
$Pr=1$
and (a)
$\theta _c=7$
, (b)
$\theta _c=14$
, for different vertical resolutions. The blue lines represent the freezing cases, and the red lines represent the melting cases. The horizontal resolution is determined by multiplying the vertical resolution by the domain aspect ratio.
Open access
