Flow regimes in a melting system composed of binary fluid: transition from penetrative convection to diffusion

Abstract

In a horizontally heated melting system, where a solid substance is subject to melting by a warmer liquid beneath, the presence of solute in the liquid introduces a complex interplay between temperature and concentration dynamics. Employing a recently developed sharp interface method (Xue et al., J. Comput. Phys., vol. 491, 2023), we conduct direct numerical simulations to investigate the transient behaviour of the system across a broad range of Rayleigh numbers and solute concentrations. Our observations reveal distinct flow regimes: at low concentrations, the system resembles a temperature-driven melting problem, characterized by vortex rolls beneath the melting interface. As the solute concentration increases, a stably stratified layer emerges beneath the interface, leading to the transition from thermal convection to penetrative convection, which resembles those flow characteristics observed in the double-diffusive convection. This shift results from the competition between the stabilizing effect induced by solute concentration gradient and the destabilizing effect caused by temperature gradient. Otherwise in the diffusion regime, characterized by very high solute concentrations, the flow becomes static due to the complete suppression of convection by the stably stratified layer. This regime further exhibits two distinct patterns: ‘melting’ and ‘dissolution’. Beyond characterizing diverse flow patterns, our study conducts a quantitative analysis, examining heat/mass transfer, melting rates, and the evolution of temperature and concentration at the interface. These insights contribute to a better understanding of the intricate interplay between temperature and solute concentration during phase change, with implications for accurately estimating melting rates in binary fluid systems.

1. Introduction

Solid–liquid phase change phenomena, encompassing processes such as melting and solidification, finds wide-ranging applications in diverse fields such as the metallurgical industry, geophysics and oceanography. In particular, when natural convection is present, the dynamics of interface patterns, such as the formation of freckles during crystal growth (Amberg & Homsy 1993), the crystallization of magma oceans (Labrosse, Hernlund & Coltice 2007), and the melting of the Antarctic Ice Sheet (Rosevear, Gayen & Galton-Fenzi 2021), are known to be highly influenced by fluid flow and the distribution of the temperature field. In this context, the composition of the system, characterized by solute concentration, plays a pivotal role, especially in the context of binary fluid systems such as binary alloys and salt water.

In the last two decades, significant attention has been directed towards the intricate coupling problem of melting and fluid flow. The pioneering investigations can be traced back to early experiments conducted by Davis, Müller & Dietsche (1984), wherein the evolution of interface patterns during melting was studied in the presence of substantial thermal convection flows in the liquid. Subsequently, the Rayleigh–Bénard convection (RBC) physical model (Ahlers, Grossmann & Lohse 2009), featuring an upper melting liquid–solid interface, emerged as a compelling framework for exploring melting problems coupled with fluid flow. The work by Vasil & Proctor (2011) contributed by deriving theoretical scaling to predict the onset of weakly nonlinear convection at small Stefan number (${St}$) number limits. Additionally, in tandem with the melting RBC, direct numerical simulations (Ulvrová et al. 2012; Rabbanipour Esfahani et al. 2018; Favier, Purseed & Duchemin 2019; Yang et al. 2023b) have provided a comprehensive map of the evolving flow patterns with increasing effective Rayleigh number. These simulations elucidate the transition path from the initiation of convection flow to merged convective rolls, culminating in turbulence. Moreover, this body of work systematically addresses the global system responses, including heat transfer, and the characteristics of the wavy interface such as height, roughness and wavelength.

It is imperative to note that the aforementioned studies focus primarily on the one-component melting RBC system, wherein the driving force is exclusively attributed to temperature-induced stratification. However, in real-world scenarios, encountering liquids composed of multi-component substances is more commonplace. In such cases, convection is not solely propelled by thermal buoyant force but is also contingent upon the inhomogeneity of the solute concentration. This multi-component nature is frequently encountered in natural settings, such as the ocean and the Earth's mantle, resulting in intriguing flow phenomena (Radko 2013). The interplay between thermal and solutal buoyancy effects gives rise to double-diffusive convection (DDC), where the temperature gradient stabilizes (destabilizes) the flow, while the concentration gradient destabilizes (stabilizes) it, manifesting in diverse parameter regimes (Yang, Verzicco & Lohse 2016; Chong et al. 2020b; Yang et al. 2022). Moreover, when a warmer liquid induces the melting of a solid, the evolution of interface patterns becomes contingent upon the distributions of both temperature and concentration fields. These two variables are strongly coupled at the interface through the Gibbs–Thomson effect, the Stefan condition, and the well-known solute-rejection relation (Davis 2001). This coupling adds a layer of complexity to the melting RBC system, as exemplified by distinctions between iceberg melting in oceans and that in pure water (Cenedese & Straneo 2023). A recent numerical study by Yang et al. (2023a) provides a preliminary insight into the interaction between the melting interface and DDC flow, specifically examining the lateral melting of an ice block in stably stratified saline water. Also recently, Rosevear et al. (2021) and Rosevear, Gayen & Galton-Fenzi (2022) demonstrated that DDC was the first-order process controlling the ocean-driven melt rate, even more important than vertical shear due to a mean flow. However, understanding this melting system remains limited due to its inherent complexity, particularly with regard to the role played by solutes in the RBC problem, which remains unclear as per current knowledge.

In the current investigation, our objective is to explore the influences of solute concentration in a horizontal melting system subjected to heating from below. The outcomes are derived through two-dimensional (2-D) numerical simulations utilizing a volume-of-fluid based sharp interface method, as detailed in our recent work (Xue et al. 2023). The adopted numerical method allows for the separate solution of temperature and solute concentration in both liquid and solid phases, by imposing the jump conditions for these variables sharply at the interface. The buoyancy-driven convection patterns are observed to depend intricately on the coupling effects between temperature and concentration. Consequently, the shape of the interface exhibits significant variations across different parameter spaces characterizing the relative strength between temperature and concentration. While acknowledging the limitations inherent in our adoption of a 2-D numerical model, it offers notable advantages. On the one hand, it enables us to explore a broader parameter space owing to the lower computational overhead compared to three-dimensional (3-D) simulations. On the other hand, previous studies have demonstrated close similarities between 2-D and 3-D thermal-driven RBC (van der Poel, Stevens & Lohse 2013; Rabbanipour Esfahani et al. 2018), as well as DDC (Chong et al. 2020a), particularly in scenarios where the Prandtl number (${Pr}$) exceeds 1, as in our present investigation. Therefore, the present 2-D numerical study, through the variation of the dimensionless strength of solutal effects represented by ${{\varLambda }}$, can already shed some light on the transition from thermal convection at small ${{\varLambda }}$ to penetrative convection at moderate ${{\varLambda }}$, and finally to pure diffusion at high ${{\varLambda }}$.

The paper is structured as follows. In § 2, we formulate the problem, delineate the range of parameters under consideration, and provide an overview of our numerical approach. Section 3 is dedicated to an in-depth discussion of the flow patterns observed in our simulations. Quantitative statistics for different types of melting dynamics are presented in § 4. Finally, concluding remarks are presented in § 5.

2. Problem statement and methodology

The physical model governing the melting of solid by the warmer binary fluid, along with the corresponding boundary conditions, is elucidated in figure 1, where figure 1(a) represents the initial condition, and figure 1(b) illustrates the melting stage. Conceptually, the melting process unfolds as follows. A binary fluid, with solute concentration denoted by $C$, is subjected to heating from below. The pure solid phase at the top, identical to the solvent in the liquid, undergoes melting due to the influence of the warmer liquid, leading to a dilution of the solute concentration in the liquid phase. The bottom wall, situated at $z = 0$, is maintained at specified temperature and concentration values, denoted as $T_{{bot}} = T_1= T_0 + \Delta T$ and $C_{{bot}} = C_0$, respectively. The top wall, positioned at $z = H$, is characterized by boundary conditions $T_{{top}} = T_0$ and $C_{{top}} = 0$. Both walls are considered to be infinitely long so that periodic conditions are applied at the left and the right boundaries.

Figure 1. Schematic illustration of the melting problem involving a binary fluid. (a) The initial configuration for the simulations. (b) The evolving stage of the system. The significance of the symbols in the images is explained in the text.

As for the initial conditions, the temperature and concentration are set to ${T = T_0}$ and $C = 0$ in the solid phase, while in the liquid phase, they are prescribed as $T(z) = T_0 + \Delta T\,(Hh_0-z)/Hh_0$ and $C = C_0$. Here, $h_0$ represents the initial dimensionless height of the liquid layer. Throughout the melting process, the temporal evolution of the interfacial topography, denoted as $\varGamma$ in the figure, is governed by the transient distribution of temperature and concentration fields in proximity to the interface. As indicated in figure 1, the temperature is continuous across the interface and its interfacial value is denoted as $T_{\mathit{\varGamma}}$. We refer to the value of the concentration from the liquid side $C^\ell$ at the interface as the interfacial concentration $C_{\mathit{\varGamma}}$. In our numerical investigations, we assume constant and identical density $\rho$ and thermal diffusivity $\alpha$ in the solid and liquid phases, with constant kinematic viscosity $\nu$ and solutal diffusivity $D$ in the liquid phase during the phase change. Specifically, the same density in two phases implies no volume expansion during the process.

We implement a rescaling strategy by normalizing the length with $H$ and time with $H^2/\alpha$. Employing the Boussinesq approximation, the dimensionless equations governing velocity, temperature and concentration in the liquid phase are expressed as

(2.1)$$\begin{gather} \partial_t\boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} ={-}\boldsymbol{\nabla} p + {Pr}\,\nabla^2\boldsymbol{u} + {Pr}\,{Ra}\,(\theta - {{\varLambda}}\phi)\boldsymbol{e}_z, \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

(2.3)$$\begin{gather}\partial_t\theta + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\theta = \nabla^2\theta, \end{gather}$$

(2.4)$$\begin{gather}\partial_t\phi + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\phi = {Le}^{{-}1}\,\nabla^2\phi, \end{gather}$$

where $\boldsymbol {u}$ represents the velocity field, $p$ is the pressure, and the unit vector in the vertical direction is denoted by $\boldsymbol {e}_z$. The term ${Pr}\,{Ra}\,(\theta - {{\varLambda }}\phi )\boldsymbol {e}_z$ represents the buoyant force along the $z$-direction generated by inhomogeneous temperature and concentration fields. Here, $\theta =(T-T_0)/\Delta T$ and $\phi =C/C_0$ denote the dimensionless temperature and concentration, respectively. It is crucial to note that temperature and solute concentration exert opposing effects on the buoyant force, a distinctive feature compared to one-component melting systems. The Prandtl number, denoted as ${Pr}$ and defined as $\nu /\alpha$, and the Lewis number, denoted as ${Le}$ and defined as $\alpha /D$, characterize the ratio of diffusive effects between momentum and temperature, and temperature and concentration, respectively. The Rayleigh number, denoted as ${Ra}$ and defined as $g\beta _T\,\Delta T\,H^3/\alpha \nu$, quantifies the strength of the thermal buoyant force. On the other hand, ${{\varLambda }} = \beta _C C_0/\beta _T\,\Delta T$ denotes the ratio of concentration to temperature effects in the induced buoyant force. Notably, the parameter ${{\varLambda }}$ is expected to significantly influence the flow pattern in the liquid, as it alters the strength of the buoyant effects. For dimensional values, $g$ represents gravity along the $z$-direction, and $\beta _{T},\beta _{C}$ are the expansion coefficients induced by temperature and concentration, respectively. It is important to highlight that in the solid phase, (2.4) becomes irrelevant since solute concentration maintains $\phi = 0$ during melting. Moreover, (2.3) transforms into a purely diffusive equation in the solid phase as the advection term is eliminated.

Compared to the one-component melting problem driven solely by the temperature, the interfacial conditions in the present binary system are more complex. The temperature at the interface is not only dependent on the melting point $\theta _M$ but also highly influenced by the local concentration of the solute. This strong coupling between interfacial temperature and solute concentration is expressed through the depression of the freezing point, as described by the Gibbs–Thomson relation (Rosevear et al. 2021; Xue et al. 2023), which can be written in dimensionless form as

(2.5)\begin{equation} \theta_\varGamma = \theta_M + m_L\phi_\varGamma. \end{equation}

In this context, the microscopic kinetic and capillary undercooling effects are considered negligible compared to the depression effect, thus are omitted from the correlation, as adopted by previous studies (Favier et al. 2019; Rosevear et al. 2021; Li, Jiao & Jia 2022). In the Gibbs–Thomson relation (2.5), $m_L$ represents the dimensionless slope of the linearized liquidus line of the binary liquid, defined as $m_L =\beta _T {\varLambda } M_L/\beta _C$, where $M_L$ is the dimensional counterpart. Moreover, the melting process of the binary fluid necessitates the simultaneous conservation of energy and mass at the interface. This implies that the heat and concentration fluxes to the interface must be balanced by the latent heat and solute fluxes during melting (Rosevear et al. 2021). Consequently, the heat flux balance at the interface is described by the dimensionless equation

(2.6)\begin{equation} v_\varGamma = {St}\,(\boldsymbol{\nabla}\theta^s - \boldsymbol{\nabla}\theta^\ell)\boldsymbol{\cdot}\boldsymbol{n}. \end{equation}

The equation is also known as the Stefan condition. Here, $v_\varGamma$ is the propagation rate of the interface, correlated to the melt rate, and $\boldsymbol {n}$ denotes the normal vector of the interface (see figure 1). The superscripts $\ell$ and $s$ denote the liquid and solid phases, respectively. The Stefan number, defined as ${St}=c^\ell _{p}\,\Delta T/\mathcal {L}$, characterizes the ratio between specific heat and latent heat, with $c_p^\ell$ representing the specific heat capacity of the liquid, and $\mathcal {L}$ the latent heat of fusion. Additionally, the concentration flux balance at the interface yields the solute-rejection relation, leading to a dimensionless formulation

(2.7)\begin{equation} \phi_\varGamma v_\varGamma = {Le}^{{-}1}\,(\boldsymbol{\nabla}\phi^s-\boldsymbol{\nabla}\phi^\ell)\boldsymbol{\cdot}\boldsymbol{n}. \end{equation}

Equations (2.5), (2.6) and (2.7) collectively describe the interfacial conditions for temperature and solute concentration, which are strongly coupled at the interface. Xue et al. (2023) emphasized that this coupling presents significant challenges in numerical simulations, rendering the solution of (2.1)–(2.7) non-trivial due to their interdependence.

For the initial conditions, the interface has dimensionless height $h_0 = 0.2$, which is uniform along the $x$-direction. The dimensionless temperature distribution is set as $\theta (z,t=0) = 1-z/h_0$ in the liquid, and $\theta (t=0) = 0$ in the solid. The solute concentration is initialized with $\phi (t=0) = 1$ in the liquid, and $\phi (t=0)=0$ in the solid, such that the solutal flux in the solid vanishes in (2.7). To satisfy the chemical equilibrium condition (2.5), the melting temperature is set to $\theta _M= -m_L$, ensuring $\theta _\varGamma (t=0) = 0$ at the interface. Due to the temporally evolving interface, we define an effective Rayleigh number ${Ra}_e(t) = {Ra}\,\bar {h}^3(t)\,(1-\bar {\theta }_\varGamma (t))$, where operator $\bar {\cdot }$ denotes the averaged value along the wavy interface. This effective Rayleigh number more accurately describes the strength of the temperature-induced buoyant force, considering the real height of the liquid layer as $H\bar {h}(t)$ rather than $H$. In the subsequent numerical investigations, we vary only ${Ra}$ and ${{\varLambda }}$ to examine the solute's influence on the melting RBC system. In the present study, physical properties are held constant in accordance with the Boussinesq approximation, and these properties closely approximate those found in ocean-driven melting systems (Rosevear et al. 2021). For example, cold saline water has molecular viscosity $\nu = 2 \times 10^{-6}\ \mathrm {m}^2\ \mathrm {s}^{-1}$, thermal diffusivity $\alpha = 1.4 \times 10^{-7}\ \mathrm {m}^2\ \mathrm {s}^{-1}$, solutal diffusivity $D = 1.3 \times 10^{-9}\ \mathrm {m}^2\ \mathrm {s}^{-1}$, thermal expansion coefficient $\beta _T = 3.8 \times 10^{-5}\,^{\circ }{\rm C}^{{-1}}$, and solutal expansion coefficient $7.8 \times 10^{-4}\ \mathrm {kg}\ \mathrm {g}^{-1}$. These dimensional properties translate to dimensionless values close to ${Pr} = 10$, ${Le} = 100$ and $m_L = -3 \times 10^{-3} {\varLambda }$, which are utilized in the current study. Furthermore, the Stefan number is maintained at ${St} = 0.1$ to ensure that the convective turnover time scale is significantly shorter than that of the melting dynamics. This approximation allows for the solid–liquid interface to be considered nearly static in the analysis of the flow behaviour. A similar Stefan number has been adopted in the previous relevant work (Favier et al. 2019; Yang et al. 2022), given that ocean-driven melting systems also exhibit low ${St}$ values. With a dimensional temperature difference $0.05\,^{\circ }\mathrm {C}$ between the top and bottom boundaries, and varying the size of the domain in $0.115 \leqslant H \leqslant 1.15\ \mathrm {m}$, the resulting Rayleigh number is in the range $10^5 \leqslant {Ra} \leqslant 10^8$, precisely the range investigated in this study. It is important to note that $Ra_e(t)$ is a dynamic parameter dependent on the time-varying height of the liquid layer $\bar {h}(t)$. The initial effective Rayleigh number falls within the range $8\times 10^3 \leqslant {Ra}_e (t = 0) \leqslant 8\times 10^6$. Furthermore, ${\varLambda }$ is controlled within the range ${10^{-4} \leqslant {\varLambda } \leqslant 10^3}$, covering a wide parameter space from a temperature-dominated regime to a solute-dominated regime.

Numerically, addressing the challenge of enforcing the jump conditions (2.5), (2.6) and (2.7) at the interface is crucial. This is achieved through the implementation of a volume-of-fluid based sharp interface method, as elaborated by Xue et al. (2023). The key aspect of this method involves separately computing the liquid and solid phases, and the conditions (2.5), (2.6) and (2.7) are imposed directly at the interface using an embedded boundary framework. The implementation utilizes the open-source solver Basilisk (Popinet 2009, 2015), and we will not specify the numerical details herein as they can be found in the above-mentioned references. A computational domain $[0,4]\times [0,1]$ with periodic boundary conditions at left wall $x=0$ and right wall $x=L=4$ is used to mimic the infinitely long domain, with spatial resolution $2048\times 512$ for situations with low-to-moderate ${Ra}$ (${Ra} \leqslant 10^7$), and $4096 \times 1024$ for higher ${Ra}$. The numerical simulation is terminated manually when the liquid fraction reaches $90\,\%$ of the domain or the interface touches the top wall. Additional details regarding the validations, and grid independence studies, can be found in Appendix A.

3. Melting dynamics and flow structures

The dynamics of the one-component melting RBC system, driven solely by temperature (i.e. at ${\varLambda } = 0$ in the current study), has been investigated systematically in existing studies (Rabbanipour Esfahani et al. 2018; Favier et al. 2019; Yang et al. 2023b). The primary features are observed when the effective Rayleigh number $Ra_e$ surpasses a well-established threshold of critical Rayleigh number ${Ra}_{{cr}}\approx 1295.78$ for the RBC with a top melting boundary (Vasil & Proctor 2011), leading to the onset of thermal convection. Concurrently, the vortex cells beneath the rising interface generate convexities and cusps along the melting interface, with the potential to merge into wider structures once their aspect ratios reach approximately $1/3$ (Favier et al. 2019).

For the melting of the binary fluid, where ${\varLambda } > 0$, it is anticipated that the flow characteristics, coupled with the involvement of interface patterns, will depend not only on the value of ${Ra}$ but also on ${\varLambda }$. Figure 2(a) illustrates the phase diagram depicting different flow regimes across a wide parameter space of $(Ra, {\varLambda })$. These regimes are identified through the present numerical simulations and include a convection regime at small ${\varLambda }$, a diffusion regime at very large ${\varLambda }$, and a layering regime in between them when $Ra$ exceeds approximately $3\times 10^6$. The corresponding flow patterns are described in figures 2(b)–2(d). These figures reveal that by enhancing the solute concentration in the liquid, the convection regime transits to the layering regime and eventually reaches the diffusion regime. Such regime transition, characterized by weakened or suppressed vortex rolls beneath the interface, originates from the absence of solute in the melted solid, leading to a reduced concentration in the upper layer of the liquid. Consequently, the concentration exhibits a negative gradient in the vertical direction, activating a stabilizing effect on the bulk flow. This is in contrast to the destabilizing effect caused by the temperature gradient (as seen in the last term of (2.1)). As a result, with increasing ${\varLambda }$, the stabilizing effect induced by solute concentration gradually dominates, transforming the purely thermal convection flow into penetrative convection flow (Ding & Wu 2021; Wang et al. 2021a). Ultimately, a stable stratified flow without convection is established. In the subsequent subsections (§§ 3.1–3.3), we will discuss each of these flow patterns in detail.

Figure 2. Flow patterns exhibited by the binary fluid melting system. (a) The distribution of these flow patterns on a phase map of $(Ra, {\varLambda })$ investigated in this study. The inset provides a magnified view of the region $1 \leqslant {\varLambda } \leqslant 10$ for improved visualization. The green dotted line at ${Ra} \approx 1.62 \times 10^5$ marks the threshold for the onset of RBC at the beginning of the simulation, estimated based on the critical Rayleigh number ${Ra}_{{cr}} \approx 1295.78$ derived by Vasil & Proctor (2011). Sketches of (b) the convection regime, (c) the layering regime, and (d) the diffusion regime are provided for clarity.

3.1. Convection regime

Let us begin by discussing the convection regime, taking the case $(Ra, {\varLambda }) = (10^6, 10^{-2})$ as an example. Figure 3(a) provides snapshots of the temperature field (left half), the concentration field (right half) and the solid–liquid interfaces (white solid lines). Clearly, a ‘cell merging’ phenomenon, similar to that observed in the melting RBC system driven solely by temperature (Rabbanipour Esfahani et al. 2018; Favier et al. 2019), is evident. The number of convective cells decreases as the interface rises, as convection spreads horizontally with the increasing liquid layer. In this case, ${\varLambda } = 10^{-2}$ is too small to induce any observable stabilizing effect on the flow. It is important to note that although a thin stably stratified layer (SSL) is initially generated beneath the interface due to the locally huge concentration gradient, it is quite vulnerable. The lower convection flow can easily penetrate this layer, leading to the formation of convexities and cusps on the interface features that are qualitatively similar to those observed at ${\varLambda } = 0$. This penetration behaviour and interface deformation are reminiscent of the findings reported by Wang et al. (2021b) for the icing problem in pure water, where density anomalies also lead to the coexistence of stable and unstable stratifications. Given that the flow structure in the convection regime has been investigated thoroughly by Favier et al. (2019), and additional numerical results for quantitative illustration are presented in Appendix A, we will not delve into further details on the flow structure in this regime.

Figure 3. Characteristics of the convection regime at $(Ra, {\varLambda }) = (10^6, 10^{-2})$. (a) Snapshots of the temperature field (left-hand image), the concentration field (right-hand image), and the solid–liquid interfaces (white solid lines) at different times: $t = 5\times 10^{-3}, 0.1, 0.3, 0.6$. (b–d) Time evolution of the horizontally averaged (b) temperature $\tilde {\theta }$, (c) concentration $\tilde {\phi }$, and (d) effective density $\tilde {\rho }$ along the vertical direction $z$. Here, $\tilde {\theta }$, $\tilde {\phi }$ and $\tilde {\rho }$ denote the averaged values at a given $z$. The red portion of a single line represents the liquid phase, while the purple portion indicates the solid phase, with the gradient depicting the presence of a wavy interface at this height. The insets of (b) and (c) present sharp temperature and concentration profiles within the slice $x=2$ at $t=0.48$, revealing interfacial values $\theta _\varGamma \approx 0$ and $\phi _\varGamma \approx 0.2$, discussed further in § 4.2.

In figure 3(b), the temporal distribution of temperature along the vertical direction is illustrated, where $\tilde {\theta }(z,t) = ({1}/{L})\int _0^L\theta (x,z,t)\,\mathrm {d}\kern0.07em x$ represents the averaged temperature at a given dimensionless height $z$ and time moment $t$. The red portion of a single line denotes the liquid phase, and the purple portion represents the solid phase, with the gradient indicating the presence of a wavy interface at this height. The dotted line represents the initial temperature distribution. Observations reveal that as RBC initiates during the melting process, the initially linearly distributed temperature transforms into a stair-like profile within the liquid. This pattern is characterized by the well-known ‘BL-shortcut-BL’ structure seen in classic RBC (Grossmann & Lohse 2000), with BL abbreviating the boundary layer. Specifically, $\tilde {\theta }$ experiences a sharp decrease in the lower and upper boundary layers, levelling off in a plateau between the two boundary layers. Notably, the value of this plateau remains nearly constant over time, estimated to be $1/2$ (indicated by the dashed horizontal line) from the numerical results. This observation aligns with findings reported by Favier et al. (2019), and we will utilize this conclusion in § 4 for quantitative analysis. Additionally, Chong et al. (2020b) observed a similar temperature profile in DDC at ${\varLambda } = 0.1$. In figure 3(c), the corresponding curves for the distribution of concentration $\tilde {\phi }$ are depicted, representing the averaged value at $z$. Similar BL-shortcut-BL profiles are observed, showing a sharp decrease at the lower and upper boundary layers in the liquid, but stabilizing at a plateau $\tilde {\phi } = 1/2$ at the developed stage. Notably, the thickness of the concentration boundary layer is smaller than that of the temperature, attributable to the solute concentration having a smaller diffusion coefficient compared to temperature, characterized by ${Le} = 100$. It is crucial to acknowledge that the $\tilde {\cdot }$ notation is employed to smooth out the discontinuity in the gradient of $\theta$ and $\phi$ at the interface. Consequently, $\tilde {\theta }$ or $\tilde {\phi }$ may not precisely capture their interfacial values at a specific $x$. To illustrate this point, we examine the vertical distribution of $\theta$ and $\phi$ at $x=2$ and time $t=0.48$, as depicted in the insets of figures 3(b) and 3(c). Here, we observe that the local interfacial temperature is approximately $\theta _\varGamma \approx 0$, while the local concentration is approximately $\phi _\varGamma \approx 0.2$. Importantly, these interfacial values remain consistent across various $x$ locations, a topic we will delve into further in § 4.2. Additionally, we estimate the averaged effective density at height $z$ using $\tilde {\rho }(z,t)={\varLambda }\,\tilde {\phi }(z,t) - \tilde {\theta }(z,t)$, and present it in figure 3(d). Here, we observe a stair-like profile that generally increases with $z$, indicating that the lighter fluid at the bottom of the liquid tends to rise upwards, triggering RBC.

However, the question arises: does the concentration field act as a passive tracer at small ${\varLambda }$? To elucidate the role that the solute plays in this scenario, figure 4(a) presents the temporal evolution of the temperature distribution located at $(x, \bar {h}(x,t)/2)$, with constant $Ra=10^6$ but ${\varLambda }$ varying between $0.1, 0.5, 1$. The white dotted line in the image indicates the first moment of observing cell merging. Clearly, a larger ${\varLambda }$ results in earlier cell merging, suggesting that the solute helps to destabilize the flow. This observation appears to be inconsistent with the previously mentioned solute-stabilizing effect. However, this apparent paradox can be understood by considering the significant difference in diffusivities between the two scalars, with ${Le} = \alpha /D = 100$. In simpler terms, when the descending flow reaches the bottom boundary layer, the melted fresh fluid carried along by the thermal plume cannot immediately adapt its concentration with the surrounding fluid due to the small diffusivity of the solute concentration. These diluted fluids are then swept horizontally, generating a fluctuating upward buoyancy force because of the lower concentration. This fluctuated motion can penetrate and agitate the thermal boundary layer at the bottom, eventually creating additional thermal plumes that further destabilize the flow. Figures 4(b) and 4(c) support this argument by focusing on the flow structures at $t = 0.228$ for ${\varLambda } = 1$. We observe thermal plumes (figure 4c) launching from positions where lower concentrations of solute are located (figure 4b), as highlighted by the closed dotted circles.

Figure 4. Influences of ${\varLambda }$ in the convection regime while maintaining the Rayleigh number at $Ra = 10^6$. (a) Temporal evolution of the temperature distribution located at $(x,\bar {h}(t)/2)$, with ${\varLambda }$ varied between $0.1$, $0.5$ and $1$. The white dotted line indicates the occurrence of the first cell merging in each case. (a) Images reveal that a larger ${\varLambda }$ leads to earlier cell merging. Contour maps of (b) solute concentration and (c) temperature are shown at $t = 0.228$ for ${\varLambda } = 1$. Comparing (b,c) indicates that thermal plumes are launched from positions where fresher plumes are located.

3.2. Diffusion regime

Now let us examine the diffusion regime. Figure 5 presents snapshots of the temperature field (left-hand images), the concentration field (right-hand images) and the solid–liquid interfaces (white solid lines), with figures 5(a,b) corresponding to the cases $(Ra, {\varLambda }) =(10^6, 5)$ and $(Ra, {\varLambda }) =(10^6, 10^3)$, respectively. In the case $(Ra, {\varLambda }) =(10^6, 5)$, it is evident that although convection rolls are still observable at $t = 10^{-2}$, their strength diminishes by $t = 5\times 10^{-2}$ and is fully suppressed by $t = 1$. As the solid melts, the fresh fluid forms an SSL that develops downwards, compressing the convection structure. However, the convection rolls are unable to penetrate the SSL. Consequently, the solute-stabilizing effect gradually dominates over the temperature-destabilizing effect, and ultimately, the fluid becomes almost static, with heat and mass transported only in a diffusive manner. In the case $(Ra, {\varLambda }) =(10^6, 10^3)$, as shown in figure 5(b) at the same time moments, we observe the convection flow transitioning to a diffusion structure as soon as the melt begins. This transition happens more rapidly than in the former case. Additionally, figure 5 reveals that the temperature exhibits stratification only in the liquid phase for ${\varLambda } = 5$, while such stratification spreads to the solid phase for ${\varLambda } = 10^3$. These distinct scenarios imply that the underlying physics of melting for ${\varLambda } = 5$ and $10^3$ are not identical, even though both belong to the diffusion regime.

Figure 5. Snapshots depicting the temperature field (left-hand images), the concentration field (right-hand images), and the positions of the solid–liquid interfaces (illustrated by white solid lines) within the diffusion regime. These snapshots correspond to scenarios at (a) $(Ra, {\varLambda }) = (10^6, 5)$ and (b) $(Ra, {\varLambda }) = (10^6, 10^3)$. The progression of time is displayed from bottom to top: $t = 10^{-2}$, $5\times 10^{-2}$ and $1$.

To elucidate this issue, let us revisit the Gibbs–Thomson correlation that couples the temperature and concentration at the interface: $\theta _\varGamma (x,t) = \theta _M + m_L\,\phi _\varGamma (x,t)$ (see (2.5)), where $m_L\,\phi _\varGamma (x,t)$ represents the concentration-induced undercooling effect. It is important to note that in our study, we have $m_L = -3\times 10^{-3}{\varLambda }$. For a small ${\varLambda }$, this implies $\theta _\varGamma (x,t) \approx \theta _M$, resulting in an adiabatic interface and an almost homogeneous temperature in the solid phase during melting. This observation aligns with what we see in figure 5(a) for ${\varLambda } = 5$. In contrast, for a much larger value of ${\varLambda }$, the undercooling effect becomes more pronounced, causing $\theta _\varGamma (x,t)$ to increase with time. This variation in $\theta _\varGamma (x,t)$ leads to a more continuous heat flux across the interface, as depicted in figure 5(b). To provide a more quantitative understanding, figure 6 illustrates the time evolution of the temperature and concentration distribution across the entire domain for different values of ${\varLambda }$ (specifically, $5$, $10^2$ and $10^3$). As we increase ${\varLambda }$, we observe the final temperature distribution aligning with our expectations: ${\varLambda } = 5$ results in $\tilde {\theta }$ distributing quasi-linearly only in the liquid, but remaining almost $0$ in the solid, while ${\varLambda } > 10^2$ corresponds to a quasi-linear $\tilde {\theta }$ in both liquid and solid. Additionally, a small ${\varLambda }$ corresponds to a more concave and nonlinear curve for $\tilde {\phi }$, while a larger ${\varLambda }$ leads to a more pronounced jump in $\tilde {\phi }$ at the interface.

Figure 6. Temporal evolution of the horizontally averaged (a,c,e) temperature and (b,d, f) concentration in the vertical direction, maintaining $Ra = 10^6$ and varying ${\varLambda }$ between $5, 10^2, 10^3$. The descriptions are consistent with those of figure 3. It is evident that at ${\varLambda } = 5$, the temperature remains nearly constant at $\tilde {\theta } = 0$ within the solid, while becoming linearly distributed across the entire domain at ${\varLambda } > 10^2$. Concerning the concentration, a higher ${\varLambda }$ leads to a more pronounced discontinuity of $\tilde {\phi }$ at the interface.

Having understood the impact of ${\varLambda }$ on temperature and concentration distributions in the diffusion regime, we revisit the Stefan condition (2.6), which governs the propagation velocity of the interface. At relatively small ${\varLambda }$ in this regime, such as ${\varLambda } = 5$, the discontinuity of the thermal flux at the interface still drives the phase change, thus it can be considered a ‘melting’ phenomenon (Woods 1992). However, in the case of a very large ${\varLambda }$, e.g. ${\varLambda } = 10^3$, the negligible jump in the thermal flux at the interface means that the phase change is no longer driven primarily by temperature but strongly depends on the solute concentration. In this case, this diffusion regime should be identified as ‘dissolution’ rather than ‘melting’, as distinguished by Woods (1992). We will distinguish these two diffusion regimes, namely melting and dissolution, in the following paragraphs because they exhibit different behaviours in mass and heat transfer.

3.3. Layering regime

Finally, let us discuss the observed layering structure at moderate values of ${\varLambda }$. The term ‘layering’ is used here to describe the stratified structures similar to those documented in the DDC literature (Rosevear et al. 2021, 2022; Yang et al. 2022, 2023a). This structure shares similarities with the flow characteristics observed in DDC scenarios. In this regime, we anticipate that the solute-stabilizing effect will prevail in the upper liquid layer, as convection is unable to fully penetrate the stably stratified layer (SSL). Conversely, the temperature-destabilizing effect will dominate in the lower liquid layer, since the SSL does not extend downwards sufficiently to completely inhibit the convection flow. As a result, an equilibrated double-layer structure may emerge, reflecting the interplay between solute and temperature effects in this intermediate regime.

Figure 7(a) illustrates a case within the layering regime at $(Ra, {\varLambda }) = (10^7, 3)$, where the coexistence of upper and lower liquid layers is evident as the solid melts. Additionally, the final distributions of temperature and concentration, as shown in figures 7(b) and 7(c), reveal a clear double-layer structure (Wang et al. 2021a; Ding & Wu 2021) within the liquid, akin to the findings reported by Rosevear et al. (2021, 2022) in their study of ocean-driven melting of horizontal ice. This structure transitions from a BL-shortcut-BL to linear diffusion from the bottom to the top of the liquid. The former represents the convective flow, while the latter characterizes the diffusive behaviour. To quantify this layering structure, we define a new dimensionless parameter, the local density ratio, given by ${\varLambda }^*(\boldsymbol {x},t) = {\varLambda }\,\partial _z \phi (\boldsymbol {x},t) / \partial _z \theta (\boldsymbol {x},t)$. We then estimate the averaged value of ${\varLambda }^*(\boldsymbol {x},t)$ at a given height, and figure 7(e) presents the spatio-temporal evolution of $\widetilde {\varLambda ^{*}}$. The figure shows that at any given time $t$, the convection and diffusion layers are distinctly separated along the height $z$, with a narrow band of high $\widetilde {\varLambda ^{*}}$ (indicated by the dark red colour) situated between them. Furthermore, the solid–liquid interface (denoted by the blue solid line) rises as the diffusion layer thickens. This ‘narrow band region’ was also identified by Rosevear et al. (2022) in their investigation of DDC in ice–ocean interactions; however, the development of the diffusion layer was not as apparent in their study. The larger melting time scale in their work resulted in a nearly stationary solid–liquid interface during the investigated period. Additionally, within the diffusion layer, $\widetilde {\varLambda ^{*}}$ is observed to decrease with height, consistent with the diminishing density gradient with height (i.e. a more flattened curve with $z$) as shown in figure 7(d). Specifically, $\widetilde {\varLambda ^{*}}$ decreases to a value less than 1 in the fluid layer adjacent to the solid–liquid interface, where $\widetilde {\varLambda ^{*}} < 1$ indicates an unstably stratified flow. This observation raises the question: can this unstable stratification in the diffusion layer trigger convection?

Figure 7. Characteristics of the layering regime with penetrative convection at $(Ra, {\varLambda }) = (10^7, 3)$. (a) Snapshots of the temperature field (left-hand image), concentration field (right-hand image) and solid–liquid interfaces (indicated by white solid lines) at times $t = 10^{-2}$, $0.1$, $0.5$ and $1.5$, from bottom to top. (b–d) The time evolutions of the vertical distribution of the averaged temperature, concentration and density at a given height, respectively. Descriptions are consistent with those provided in figure 3. (e) The time evolution of the averaged local density ratio along the height, $\widetilde {{\varLambda }^*}(z,t)$, in the $(t,z)$ plot, with the blue solid line denoting the position of the solid–liquid interface.

We investigate this by increasing ${Ra}$ to a higher value, such as $({Ra}, {\varLambda }) = (10^8, 3)$, as shown in figure 8. Figure 8(a) illustrates a convection–convection double-layer structure in the liquid. Here, the unstable stratification below the interface destabilizes the diffusion layer, leading to its transition into the upper convection layer. This double-layer structure resembles the well-known ‘thermohaline staircase’ observed in the DDC literature (Rosevear et al. 2021; Yang et al. 2020, 2022). It is characterized by the double staircases in the profiles of $\tilde {\theta }$ and $\tilde {\phi }$, as depicted in figures 8(b) and 8(c). As presented in figure 8(e), the two layers are separated by a narrow band of high $\widetilde {{\varLambda }^*}$, whose strength decreases along the melting process. This implies a further destabilization of this double-layer structure. Although this transition is not captured in this case due to the top boundary, we have further simulated a case $({Ra}, {\varLambda }) = (10^7, 2.5)$ with an enlarged domain in Appendix B, which allows us to capture the entire transition from a convection–diffusion layering regime to a convection regime. Additionally, figure 2(a) demonstrates that for $Ra < 3 \times 10^6$, the layering regime does not manifest in the simulations. These observations strongly suggest that the layering regime is transient, ultimately evolving into either a convection-dominant regime or a diffusion-dominant regime as melting progresses.

Figure 8. Characteristics of the layering regime with convection–convection double-layer structure at $(Ra, {\varLambda }) = (10^8, 3)$, with descriptions identical to those in figure 7, except that the time moments of the snapshots in (a) are $t=10^{-2}$, 0.1, 0.3 and $0.5$ from bottom to top.

We now revisit figure 2(a), focusing on the region of relatively low Rayleigh number ($Ra$), specifically below the green dotted line that denotes the threshold for the onset of RBC at the beginning of the simulation. In this regime, both the temperature and concentration fields initially exhibit purely diffusive transport. The onset of convection, as the solid melts, is influenced by local stratification. According to Radko (2013), a critical value ${\varLambda }_{{cr}} = ({Pr}+1)/({Pr}+{Le}^{-1}) \approx 1.1$ is predicted for the onset of thermohaline convection through linear stability analysis (LSA). This critical value is slightly higher than our numerical result ${\varLambda }_{{cr}} \approx 1.0$. This discrepancy is attributed to the fact that LSA assumes a linear concentration profile, whereas our study observes a concave profile due to its lower diffusivity. This nonlinearity leads to a larger local density ratio, making convection more difficult to initiate. In contrast, examining results in the higher ${Ra}$ region, above the green dotted line, reveals that an increased ${Ra}$ promotes the convection regime. This observation indicates that at higher Rayleigh numbers, the threshold ${\varLambda }_{{cr}}$ for transitioning from the convection regime to the layering regime increases. The observed increase in ${Lambda}_{{cr}}$ can be attributed to the dynamics of penetrative convection, where a convection-dominated fluid layer extends into an initially stably stratified layer (an SSL). As demonstrated by Wang et al. (2021a,b), a higher Rayleigh number enhances the heat transfer within the convection layer, necessitating a reduction in the thickness of the SSL to balance the increased heat flux from the convection layer. Consequently, with higher ${Ra}$, the SSL becomes thinner and more susceptible to penetration. Under these conditions, a higher ${\varLambda }$ is required to provide a stronger stabilizing effect, which is essential for maintaining the equilibrium structure of the layering regime at elevated ${Ra}$ values.

It is also noteworthy to discuss the unusual distribution of $\tilde {\rho }$ within the convection layers of the layering regime. Contrary to the ascending density profile observed in figure 3(d) for $\varLambda = 10^{-2}$, the density profiles presented in figures 7(d) and 8(d) exhibit a descending, stair-like pattern. This suggests that lighter fluid is situated above heavier fluid, which would typically indicate a stably stratified flow rather than an unstable convective flow. To resolve this inconsistency, the insets in figures 7(d) and 8(d) display magnified views of the density distribution within the bulk liquid, revealing a gravitationally unstable stratified region characterized by increasing density with height. At $t=0$, the convective motions are initiated by this initially unstable stratification, as indicated by the dotted line in figures 7(d) and 8(d). Over time, the solutal boundary layer evolves and forms a thick, stably stratified boundary layer, while the bulk liquid retains its unstable stratification. This phenomenon is consistent with observations reported by Chong et al. (2020b), who describe it as a subcritical behaviour in DDC.

4. Quantitative analysis

After examining various flow patterns resulting from the manipulation of ${\varLambda }$ and $Ra$, we delve into a more quantitative analysis of the flow characteristics and interfacial behaviours. This encompasses a detailed examination of heat and mass transfer, the interfacial values of temperature and concentration, and the evolution of the interface height.

4.1. Heat and mass transfer

To assess heat transfer in a melting system, a crucial dimensionless parameter is the thermal Nusselt number, denoted as ${Nu}_T$, which characterizes the ratio of heat transfer caused by convection to conduction. It is important to note a distinction in ${Nu}_T$ between the RBC system with and without melting. In classical RBC, the entire system is in energetic equilibrium, allowing ${Nu}_T$ to be calculated from any horizontal plane. In the presence of melting, however, the injected energy must include latent heat due to melting, resulting in spatial variations in ${Nu}_T$. Here, we define ${Nu}_T$ using the injected flux, expressed as

(4.1)\begin{equation} {Nu}_T(t) = \frac{Q_C(t)}{Q_D(t)} = \frac{\displaystyle\int\nolimits_{0}^L \left.-\dfrac{\partial\theta}{\partial z}\right|_{z=0}(x, t)\,\mathrm{d}\kern0.07em x}{L\,\dfrac{1-\bar{\theta}_\varGamma(t)}{\bar{h}(t)} \left(1+\dfrac{2\,{\Delta A}(t)}{\rm \pi}\right)}, \end{equation}

where $Q_C$ and $Q_D$ are the injected thermal flux realized by the convection and by the pure diffusion, respectively. We recall that $x=L=4$ is the position of the right-hand wall where periodic boundary condition is imposed. Besides, the term $2\,{\Delta A(t)}/{\rm \pi}$ in (4.1) represents a second-order correction considering topography, with $\Delta A = A - 1$, and $A$ the area of the solid–liquid interface per unit length (Favier et al. 2019). Figure 9(a) illustrates the evolution of ${Nu}_T(t)$ with respect to $Ra_e(t)$, showing a growth with melting, as explained earlier. Notably, for a melting system driven solely by temperature (${\varLambda } = 0$), the thermal Nusselt number follows a scaling ${Nu}_T \sim Ra_e^{1/4}$, as confirmed by plotting the result of $(Pr, St) = (10, 0.1)$ from Rabbanipour Esfahani et al. (2018). In our study, we find that for cases with ${\varLambda } \leqslant 1$ in the convection regime, all curves ultimately collapse onto the correlation ${{Nu}_T} \sim {Ra}_e^{1/4}$ after initial jumps due to the onset of convection. Furthermore, they share identical prefactors $0.25$, indicating that the solute concentration, despite varying with ${\varLambda }$, has minimal influence on the heat transfer of the system. However, when ${\varLambda }$ approaches the borderline between convection and layering regimes (e.g. ${\varLambda }=2$ at ${Ra}=10^7$ indicating weak convection), the stabilizing solutal effect becomes non-negligible for heat transfer. Specifically, we observe ${Nu}_T$ initially dropping but eventually converging to ${Ra}_e^{1/4}$ after a sharp rise. This initial drop is attributed to the formation of the stabilizing effect from the SSL temporarily separating the convective liquid layer from the solid–liquid interface. However, as the interface height rises during melting, $Ra_e$ increases, allowing convection to penetrate the SSL, deform the interface, and sharply enhance heat transfer. Similarly, the ${Nu}_T$ profile corresponding to the layering regime (e.g. ${\varLambda } = 3$ at $Ra= 10^8$) exhibits a variation trend similar to that of the weak convection regime, but with exponent values smaller than $1/4$. Moreover, for the diffusion regime with large ${\varLambda }$, the convection flow is entirely eliminated by the growth of SSL. Consequently, ${Nu}_T$ initially jumps to a large value due to the onset of convection, but ultimately converges to 1 in a short time period.

Figure 9. The temporal evolution of (a) the thermal Nusselt number ${Nu}_T(t)$ and (b) the solutal Nusselt number ${Nu}_C(t)$, plotted against ${Ra}_e(t)$ for varying $Ra$ (${Ra}=10^6,10^7,10^8$) and ${\varLambda }$ ($10^{-2} \leqslant {\varLambda } \leqslant 10^3$). In (b), the dash-dotted line represents the predicted value of ${Nu}_C\approx 2.49$ for the ‘melting’-dominated diffusion regime (${\varLambda } \leqslant 10$) obtained through scaling analysis. For better readability, the different groups of profiles with ${Ra}=10^6$, $10^7$ and $10^8$ are distinguished by different transparency. Also, the range of ${Ra}_e$ for different ${Ra}$ is highlighted in both plots.

To evaluate mass transfer in this system, we introduce the solutal Nusselt number, denoted as ${Nu}_C$. This parameter is defined in a manner analogous to the thermal Nusselt number ${Nu}_T$ and is expressed as

(4.2)\begin{equation} {Nu}_C(t) = \frac{\displaystyle\int\nolimits_{0}^L \left.-\dfrac{\partial\phi}{\partial z}\right|_{z=0}(x, t)\, \mathrm{d}\kern0.07em x}{L\,\dfrac{1-\bar{\phi}_\varGamma(t)}{\bar{h}(t)}\left(1+\dfrac{2\,{\Delta A}(t)}{\rm \pi}\right)}, \end{equation}

which represents the ratio between the mass flux induced by convective processes and that resulting from pure diffusion. Figure 9(b) presents the results, noting that since the binary fluid has an initially homogeneous concentration $\phi = 1$, independent of ${\varLambda }$, all ${Nu}_C$ values start from 0. In the convection and layering regimes, the evolution of ${Nu}_C$ closely mirrors that of ${Nu}_T$ because both the injected thermal and mass fluxes are influenced by the same convective structures. In particular, the solutal Nusselt number converges to ${Nu}_C = 1.4\,{Ra}_e^{1/4}$ at ${\varLambda } \leqslant 1$, while it is ${Nu}_T = 0.25\,{Ra}_e^{1/4}$ for the thermal Nusselt number. It is crucial to understand the ratio of the two prefactors, i.e. why ${Nu}_C(t)/{Nu}_T(t) \approx 1.4/0.25 = 5.6$. A scaling analysis based on the energy/mass budget helps to clarify this ratio.

Regarding heat transfer, the injected thermal flux serves two primary functions. First, it maintains the liquid at an average temperature, denoted by $\langle \theta \rangle _\ell \approx 1/2$, as illustrated in figure 3(b), where $\langle \cdot \rangle _\ell$ represents the average over the liquid phase. Second, it provides the necessary heat flux to facilitate the melting of the solid. These contributions are quantified by $({L}/{2}) ({\mathrm {d}\bar {h}(t)}/{\mathrm {d}t})$ and $\int _{s\unicode{x2013}l\ interface} [(- \boldsymbol {\nabla }\theta ^\ell )\boldsymbol {\cdot }\boldsymbol {n}](s,t)\,\mathrm {d}s$, respectively, with $\int _{s\unicode{x2013}l\ interface}$ indicating the line integral of the heat flux along the solid–liquid interface. Consequently, the energy budget can be expressed as

(4.3) \begin{equation} \int_{0}^L \left.-\frac{\partial\theta}{\partial z}\right|_{z=0}(x, t)\,\mathrm{d}\kern0.07em x \approx \frac{L}2\,\frac{\mathrm{d}\bar{h}(t)}{\mathrm{d}t} + \int_{s\unicode{x2013}l\ interface} [(- \boldsymbol{\nabla}\theta^\ell)\boldsymbol{\cdot}\boldsymbol{n}](s,t)\,\mathrm{d}s. \end{equation}

Substituting (4.1) and (2.6) into this equation further reformulates it to

(4.4)\begin{equation} \frac{{Nu}_T(t)\,L}{\bar{h}(t)} \approx \frac{L}2\,\frac{\mathrm{d}\bar{h}(t)}{\mathrm{d}t} + {St}^{{-}1}\int_{s\unicode{x2013}l\ interface} v_\varGamma(s,t) \,\mathrm{d}s, \end{equation}

where several approximations are adopted: (1) the change in the interface area is negligible, i.e. $\Delta A \ll 1$; (2) the averaged interfacial temperature in the convection regime is nearly zero, i.e. $\bar {\theta }_\varGamma \approx 0$, as discussed in § 4.2; (3) the heat flux contribution from the solid phase is insignificant, i.e. $(\boldsymbol {\nabla }\theta ^s)\boldsymbol {\cdot }\boldsymbol {n} \approx 0$, as identified in figure 3(b). Additionally, the integral of the interface propagation rate, which appears as the second term on the right-hand side of (4.4), can be estimated based on the average height as

(4.5)\begin{equation} \int_{s\unicode{x2013}l\ interface}v_\varGamma(s,t)\,\mathrm{d}s = L\,\frac{\mathrm{d}\bar{h}(t)}{\mathrm{d}t}. \end{equation}

Finally, (4.3) becomes

(4.6)\begin{equation} \frac{{{Nu}_T}(t)}{\bar{h}(t)} \approx \frac{\mathrm{d} \bar{h}(t)}{\mathrm{d}t}\left(\frac12 + {St}^{{-}1}\right), \end{equation}

which describes the dynamic coupling between the heat transfer and the evolution of the solid–liquid interface.

Regarding the mass budget, the injected mass flux serves to maintain the liquid at an average concentration $\langle \phi \rangle _\ell \approx 1/2$, as shown in figure 3(c). The mass budget can be expressed as

(4.7)\begin{equation} {Le}^{{-}1}\int_{0}^L \left.-\frac{\partial\phi}{\partial z}\right|_{z=0}(x, t)\,\mathrm{d}x\approx \frac{L}2\,\frac{\mathrm{d}\bar{h}(t)}{\mathrm{d}t}. \end{equation}

Similarly, by substituting (4.2) and applying the approximation $\Delta A \ll 1$, the injected mass flux can be reformulated as

(4.8)\begin{equation} {Le}^{{-}1}\,\frac{{Nu}_C(t)\,(1-\phi_\varGamma(t))}{\bar{h}(t)} \approx \frac12\, \frac{\mathrm{d}\bar{h}(t)}{\mathrm{d}t}, \end{equation}

which reflects the interplay between mass transfer and the evolving interface. Combining the energy and mass conservation laws given by (4.6) and (4.8), we derive the expression

(4.9)\begin{equation} \frac{{Nu}_C(t)}{{{Nu}_T}(t)} \approx \frac{{Le}}{(1+2\,{St}^{{-}1})(1- \bar{\phi}_\varGamma(t))} \approx 6, \end{equation}

where $\bar {\phi }_\varGamma (t)$ is known to be approximately 0.2 when ${\varLambda } \leqslant 1$ (see figure 3 and the corresponding text), and this relationship will be discussed further in § 4.2. The predicted value ${Nu}_C(t)/{Nu}_T(t) \approx 6$ closely matches the numerical observation $5.6$. This analysis provides a quantitative understanding of the ratio between solutal and thermal Nusselt numbers in the context of mass and energy conservation.

In the diffusion regime, which is categorized into ‘melting’ or ‘dissolution’ patterns, we can provide further insights. For the dissolution cases (${\varLambda } \geqslant 10^3$), where the $\phi$ profile is almost linearly distributed in the liquid (see figures 6b,d, f), we expect the solutal Nusselt number to approach 1, similar to the thermal Nusselt number in the diffusion regime. However, for the melting cases, the situation becomes more complex as the concave $\phi$ profile in the liquid exhibits more nonlinearity (see figures 6b,d, f). In this scenario, we establish a theoretical argument by approximating the entire process as a fully melting one-dimensional Stefan problem, enabling us to obtain self-similar solutions. We assume that $\theta _\varGamma$ and $\phi _\varGamma$ are maintained at 0 during melting, as illustrated in figure 6. The temperature and concentration fields in the liquid can then be formulated as $\theta (z,t)=1-\mathrm {erf}(z/2\sqrt {t})/\mathrm {erf}(\xi )$ and $\phi (z,t)=1-\mathrm {erf}(z\sqrt {{Le}}/2\sqrt {t})/\mathrm {erf}(\xi \sqrt {{Le}})$. The interface evolves as $h(t)=2\xi \sqrt {t}$, with $\xi$ being obtained from the heat and mass balance at the interface, given that ${St} = \xi \,\mathrm {erf}(\xi )\exp (\xi ^2)\sqrt {{\rm \pi} }$. Note that these solutions have been discussed in our previous study (Xue et al. 2023). Based on the presented arguments, the solutal Nusselt number can be predicted as ${Nu}_C = -\partial _z\phi |_{z=0}\,h(t)\approx 2.49$. Consequently, for all cases in the diffusion regime, we would expect the stabilized ${Nu}_C$ to fall in a region $1 < {Nu}_C < 2.49$, with the two extremes being labelled by the dashed lines and dash-dotted lines in figure 9(b). This prediction aligns with our numerical results at ${\varLambda } \geqslant 5$ in the same plot.

4.2. Interfacial temperature and concentration

In this subsection, we investigate the temporal evolution of the averaged temperature and concentration at the interface, denoted by $\bar {\theta }_\varGamma (t)$ and $\bar {\phi }_\varGamma (t)$, respectively. We emphasize their strong coupling through the Gibbs–Thomson correlation (2.5) that $\bar {\theta }_\varGamma (t)=\theta _M + m_L\,\bar {\phi }_\varGamma (t)$, where $m_L = -3\times 10^{-3}{\varLambda }$.

For small to moderate ${\varLambda }$, where ${\varLambda } \leqslant 10$, the average interfacial temperature $\bar {\theta }_\varGamma (t)$ remains at $0$, largely due to the negligible solutal undercooling effect. However, as is evident from figure 3, the interfacial concentration appears to stabilize around $\phi _\varGamma \approx 0.2$ for low ${\varLambda }$, warranting further investigation. Figures 10(a)–10(c) present the time histories of $\bar {\phi }_\varGamma (t)$ at different ${\varLambda }$ (${\varLambda } \leqslant 10$), with ${Ra}=10^6$, $10^7$ and $10^8$ depicted in separate plots. Additionally, figure 10(d) showcases their $\bar {\phi }_\varGamma$ values at the final moment. As melting commences, $\bar {\phi }_\varGamma (t)$ decreases from $1$, and notably, for cases where ${\varLambda } < 1$, $\bar {\phi }_\varGamma (t)$ almost stabilizes at a constant value $0.2$, independent of ${\varLambda }$. Even at ${Ra} = 10^7$ and $10^8$, which exhibit secondary convection in the upper liquid layer, $\bar {\phi }_\varGamma (t)$ converges to approximately $0.2$. This consistency suggests that the ultimate $\bar {\phi }_\varGamma$ significantly depends on the convective structure near the solutal boundary layer beneath the interface, irrespective of whether the liquid comprises one layer or double layers.

Figure 10. Time evolution of $\bar {\phi }_\varGamma (t)$ within the range of small to moderate ${\varLambda }$ (${\varLambda } \leqslant 10$) depicted at (a) ${Ra=10^6}$, (b) ${Ra}=10^7$, and (c) ${Ra}=10^8$, respectively. The corresponding ${\varLambda }$ values are represented by different line colours, consistent with those in figure 9. (d) The ultimate value of $\bar {\phi }_\varGamma$ as a function of ${\varLambda }$ for different ${Ra}$. (e) A schematic representation of the one-dimensional BL-shortcut-BL structure for the concentration $\phi$ and temperature $\theta$ profiles. The upper and lower solutal (thermal) boundary layers are symbolized by $\delta _{\phi,u}$ and $\delta _{\phi,\ell }$ ($\delta _{\theta,u}$ and $\delta _{\theta,\ell }$), respectively.

To elucidate why $\bar {\phi }_\varGamma (t)$ approaches $0.2$, we establish a theoretical model based on two assumptions. First, inspired by the profiles within the slice $x=2$ (insets of figures 3b,c), we approximate the concentration and temperature profiles in the liquid using a one-dimensional BL-shortcut-BL structure (figure 10e) observed in the convection regime. We assume $\phi$ to be linearly distributed in the boundary layer, with value $1/2$ in the bulk flow. The thicknesses of the upper and lower solutal boundary layers are given by

(4.10a,b)\begin{align} \delta_{\phi,u}(t) = \left.\left(\frac{1}{2} - \phi_\varGamma(t)\right)\right/\left(\left.-\frac{\partial\phi}{\partial z}\right|_{z=\bar{h}}(t)\right) \quad \mathrm{and} \quad \delta_{\phi,\ell}(t) = \left.\frac{1}{2}\right/\left(\left.-\frac{\partial\phi}{\partial z}\right|_{z=0}(t)\right). \end{align}

The ratio between the thicknesses of the lower and upper solutal boundary layers is given by

(4.11)\begin{equation} \frac{\delta_{\phi,\ell}(t)}{\delta_{\phi,u}(t) } = \frac{\left.\dfrac{\partial\phi}{\partial z}\right|_{z=\bar{h}}(t)}{\left(1-2\phi_\varGamma(t)\right)\left.\dfrac{\partial\phi}{\partial z}\right|_{z=0}(t)}, \end{equation}

whereas the upper mass flux can be estimated using the solute-rejection relation (2.6), with the approximation ${\partial \phi }/{\partial z}|_{z=\bar {h}}(t) \approx -{Le}\,\phi _\varGamma (t)\,({\mathrm {d}\bar {h}(t)}/{\mathrm {d}t})$, and the lower mass flux is estimated from the mass budget (4.7), with ${\partial \phi }/{\partial z}|_{z=0}(t) \approx -({{Le}}/{2})({\mathrm {d}\bar {h}(t)}/{\mathrm {d}t})$. Substituting these approximations into (4.11) yields

(4.12)\begin{equation} \frac{\delta_{\phi,\ell}(t)}{\delta_{\phi,u}(t)} \approx \frac{2\phi_\varGamma(t)}{1-2\phi_\varGamma(t)}. \end{equation}

Following a method analogous to (4.10a,b), the thicknesses of the upper and lower temperature boundary layers are estimated as

(4.13a,b)\begin{equation} \delta_{\theta,u}(t) = \left.\frac{1}{2}\right/\left(\left.-\frac{\partial\theta}{\partial z}\right|_{z=\bar{h}}(t)\right) \quad\mathrm{and}\quad \delta_{\theta,\ell}(t) = \left.\frac{1}{2}\right/\left(\left.-\frac{\partial\theta}{\partial z}\right|_{z=0}(t)\right), \end{equation}

thus the thickness ratio is

(4.14)\begin{equation} \frac{\delta_{\theta,\ell}(t)}{\delta_{\theta,u}(t)} = \left.\left(\left.-\frac{\partial{\theta}}{\partial z}\right|_{z=\bar{h}}(t)\right)\right/ \left(\left.-\frac{\partial{\theta}}{\partial z}\right|_{z=0}(t)\right). \end{equation}

The upper heat flux is related to the Stefan condition (2.6), so we have $({\partial {\theta }}/{\partial z})|_{z=\bar {h}}(t) \approx {St}^{-1}\, ({\mathrm {d}\bar {h}(t)}/{\mathrm {d}t})$, by neglecting the heat flux contribution from the solid side. Additionally, the lower heat flux is determined from the energy budget (4.3) and (4.6), which gives $-({\partial {\theta }}/{\partial z})|_{z=0}(t) \approx ({\mathrm {d} \bar {h}(t)}/{\mathrm {d}t})(\frac 12 + {St}^{-1})$. Combining these two correlations with (4.14) leads to

(4.15)\begin{equation} \frac{\delta_{\theta,\ell}(t)}{\delta_{\theta,u}(t)} \approx \frac{2}{{St}+2}. \end{equation}

Additionally, it is reasonable to assume that the thickness ratio between the lower and upper boundary layers is the same for both the concentration and temperature fields. Thus we approximate that

(4.16)\begin{equation} \frac{\delta_{\phi,\ell}(t)}{\delta_{\phi,u}(t)} \approx \frac{\delta_{\theta,\ell}(t)}{\delta_{\theta,u}(t)}, \end{equation}

which can be further simplified by considering (4.12) and (4.15), leading to the correlation

(4.17)\begin{equation} \phi_\varGamma \approx \frac{1}{{St}+4} \approx 0.24. \end{equation}

This value closely matches our numerical solution. For other cases with ${\varLambda } \leqslant 10$, including the ‘melting’-dominated diffusion regime and the layering regime without secondary convection, the concentration gradually drops and converges to $0$ due to the slower diffusion of concentration compared to temperature in the diffusion layer beneath the melting interface.

For cases where ${\varLambda } \geqslant 10^2$, dominated by the ‘dissolution’ process, an intriguing interplay between $\bar {\theta }_\varGamma (t)$ and $\bar {\phi }_\varGamma (t)$ emerges. The solid lines in figure 11 illustrate the evolution of $\bar {\phi }_\varGamma$ and $\bar {\theta }_\varGamma$ with respect to the average interface height $\bar {h}$ for ${\varLambda } = 10^2$ and $10^3$. Due to a significant solutal undercooling effect, the decrease in $\bar {\phi }_\varGamma$ leads to a noticeable increase in $\bar {\theta }_\varGamma$, ultimately converging to the guide lines in the figure during the fully developed stage. Understanding this convergence is linked to the characteristics discussed in figure 6. These characteristics reveal a quasi-linear temperature profile across the entire domain in the fully developed stage, allowing the interfacial temperature $\bar {\theta }_\varGamma$ to be approximated as

(4.18)\begin{equation} \bar{\theta}_\varGamma(t) \approx 1 - \bar{h}(t). \end{equation}

Considering the solutal undercooling effect (2.5), the interfacial concentration $\bar {\phi }_\varGamma$ can be approximated as

(4.19)\begin{equation} \bar{\phi}_\varGamma(t) \approx \frac{1-\bar{h}(t)}{m_L} + 1. \end{equation}

Expressions (4.18) and (4.19) are presented as dotted lines in figure 11 for both cases, and demonstrate excellent agreement with our numerical solution in the developed stage.

Figure 11. Evolution of the averaged interfacial values of temperature $\bar {\theta }_\varGamma$ and concentration $\bar {\phi }_\varGamma$ in the ‘dissolution’-dominated regime, with the Rayleigh number maintained at ${Ra}=10^6$, for (a) ${\varLambda } = 10^2$, and (b) ${\varLambda } = 10^3$. Solid lines represent numerical results, while the dotted and dash-dotted lines depict the analytical solutions of (4.18) and (4.19), respectively.

4.3. Interface evolution

We then shift our focus to the time evolution of the interface height $\bar {h}(t)$, analysing it through an energy argument, a method employed previously by others (Rabbanipour Esfahani et al. 2018; Favier et al. 2019). In both the convection regime and the ‘melting’-dominated diffusion regime, we make two key approximations. First, we assume that the average temperature in the liquid is $\langle \theta \rangle _\ell = 1/2$. Second, we assume that $\theta _\varGamma (x,t)$ remains at 0. We define a universal correlation ${Nu}_T(t)=\gamma _1\,{Ra}_e^{\gamma _2}(t)$, with $(\gamma _1,\gamma _2) = (0.25,1/4)$ for the convection regime, and $(\gamma _1, \gamma _2) = (1, 0)$ for the ‘melting’-dominated diffusion regime, from the observation in figure 9(a). Substituting this correlation into (4.6) yields an ordinary differential equation (ODE)

(4.20)\begin{equation} \gamma_1\,{Ra}^{\gamma_2}\,\bar{h}^{3\gamma_2-1}(t)\approx\frac{\mathrm{d}\bar{h}(t) }{\mathrm{d}t}\left(\frac12+{St}^{{-}1}\right), \end{equation}

with the relation ${Ra}_e(t) = {Ra}\,\bar {h}^3(t)$. We then obtain the solution of (4.20), which reads

(4.21)\begin{equation} \bar{h}(t) \approx \left[ \frac{(2-3\gamma_2)\gamma_1\,{Ra}^{\gamma_2}}{\dfrac12+{St}^{{-}1}}\,t + h_0^{2-3\gamma_2}\right]^{1/({2-3\gamma_2})}. \end{equation}

For $Ra = 10^6$, the predicted behaviours for $\bar {h}(t)$ based on (4.21) are

(4.22)\begin{equation} \bar{h}(t) \approx \left\{\begin{array}{@{}ll@{}} (0.941t+0.134)^{4/5} & \text{for convection regime}, \\ (0.190t+0.040)^{1/2} & \text{for `melting'-dominated diffusion regime}. \end{array} \right. \end{equation}

Consequently, the averaged interface height propagates exponentially, with scaling exponents $4/5$ and $1/2$ corresponding to the convection regime and the ‘melting’-dominated diffusion regime, respectively. Accordingly, the melting rate of the solid follows the scaling exponents $-1/5$ and $-1/2$ in these two regimes. In figure 12(a), a comparison is presented between the theoretical predictions (solid black line and dashed black line for the two options of (4.22)) and the numerical solution in these two regimes, showing excellent agreement. Additionally, the scaling exponents $1/2$ and $4/5$ are plotted as grey dash-dotted lines, providing an accurate prediction of the development trend of the corresponding numerical solutions in the two regimes.

Figure 12. Characteristics of the interface height during melting. (a) Averaged interface height at ${Ra}=10^6$. (b) Roughness of the interface height at ${Ra}=10^6$, $10^7$ and $10^8$. The colours of the lines correspond to different values of ${\varLambda }$ and are consistent with those described in figure 9. In (a), the black solid line represents the theoretical predictions from (4.22) for the convection regime, the black dashed line from (4.22) for the ‘melting’-dominated diffusion regime, and the black dotted line from (4.25) for the diffusion-dominated regime. The inset in (a) presents the numerical solution for the time evolution of the average height $\bar {h}$ over an extended period for ${\varLambda } = 10^3$, illustrating the very slow dissolution process, which aligns well with the prediction from (4.25) for the diffusion-dominated regime.

To estimate $\bar {h}(t)$ in the ‘dissolution’-dominated diffusion regime (e.g. at ${\varLambda } = 10^3$), the assumption $\theta _\varGamma = 0$ is no longer valid, as seen in figures 6(a,c,e). However, it is still possible to approximate it with the asymptotic correlation between the average interfacial concentration $\bar {\phi }_\varGamma$ and the average interface height $\bar {h}$ of (4.19):

(4.23)\begin{equation} \bar{\phi}_\varGamma(t) \approx \frac{1-\bar{h}(t)}{m_L} + 1 = \frac{1}{3}\,\bar{h}(t) + \frac{2}{3}. \end{equation}

Recalling that $m_L = -3 \times 10^{-3}{\varLambda } = -3$, the ‘dissolution’ dominated by the solute $\phi$ allows us to approximate the $\phi$ profile in the liquid with the one-dimensional linear function

(4.24)\begin{equation} \phi(z,t) \approx \frac{\bar{\phi}_\varGamma(t) - 1}{\bar{h}(t)}\,z + 1. \end{equation}

Combining (4.23), the vertical gradient of $\phi$ from (4.24), and the solute-rejection relation (2.6), we obtain the ODE

(4.25)\begin{equation} {Le}^{{-}1}\,\frac{1 - \left(\dfrac{1}{3}\bar{h}(t) + \dfrac{2}{3}\right)}{\bar{h}(t)} \approx \left(\frac{1}{3}\,\bar{h}(t) + \frac{2}{3}\right)\frac{\mathrm{d}\bar{h}(t)}{\mathrm{d}t}. \end{equation}

By solving (4.25), the analytical solution for the interface height $\bar {h}(t)$ at ${\varLambda } = 10^3$ is plotted in figure 12(a), showing good agreement with the numerical solution, except for a slight underestimation at the initial stage. The deviation arises from the faster melting by the initial triggered convection, which is not considered by the theoretical estimation.

Another relevant characteristic that requires investigation is the roughness of the interface, which is characterized as

(4.26)\begin{equation} h_{{rms}}(t)=\sqrt{\overline{(h(x,t)-\bar{h}(t))^2}}\,{Ra}^{1/3} = \sqrt{\overline{h'(x,t)^2}}\,{Ra}^{1/3}, \end{equation}

where ${Ra}^{1/3}$ is used for normalization at different ${Ra}$ (Yang et al. 2023b). The numerical results are presented in figure 12(b). We observe that in all cases, $h_{{rms}}$ remains at $0$ until the onset of convection deforms the melting interface, as described previously. In the convection regime corresponding to very small ${\varLambda }$, the developing trend of the curves follows that of $h_{{rms}} \sim {Ra}_e^{1/3}$, consistent with what was derived by Yang et al. (2023b) for the pure temperature-driven melting problem. For the diffusion and layering regimes characterized by promoting ${\varLambda }$, the interface roughness is observed to reduce in the picture. This reduction is attributed to the weakening or disappearance of convection rolls in these two regimes. Additionally, even if there are cusps and convexities along the interface, the cusp regions melt more quickly since the local temperature has a negative gradient there, resulting in a decreasing $h_{{rms}}$.

5. Conclusions and outlook

In this study, we explore the impact of solute presence in a horizontally heated melting system through numerical simulations employing a recently developed sharp interface method. The system comprises a binary fluid situated beneath a pure solid substance, with the latter serving as the solvent for the liquid. As the warmer fluid beneath melts the solid from the bottom, the liquid experiences dilution from the top. The interplay between temperature and solute concentration leads to contrasting effects on buoyant forces, owing to the distinct variations in density. While temperature induces well-known Rayleigh–Bénard convection (RBC), the solute concentration acts as a stabilizing factor. Our primary objective is to delve into their antagonistic dynamics across a broad parameter space defined by $(Ra, {\varLambda })$.

Our numerical findings reveal three distinct flow regimes. The first, characterized by small ${\varLambda }$, corresponds to the convection regime, where vortex rolls form beneath the melting interface due to the relatively weak stabilizing effect induced by the solute. In addition, contrary to the initial expectation that the solute would act as a passive tracer or exhibit a stabilizing effect in this low ${\varLambda }$ scenario, it instead has a destabilizing effect on the triggered convection. This is due to the fact that the lighter fluid carries a lower solute concentration, which agitates the lower thermal boundary layer and generates a fluctuating upward buoyancy force. As ${\varLambda }$ increases, the system undergoes a transition from thermal convection to penetrative convection, giving rise to the layering regime. This regime is characterized by moderate solutal effects resulting from the melting solid, which leads to the formation of a stably stratified layer that inhibits the full dominance of convection in the liquid. In this regime, a double-layer flow structure emerges: the lower layer remains in the convection regime, while the upper layer may be either another convection regime or a distinct diffusion regime, depending on the value of ${\varLambda }$. This layering structure exhibits behaviours very similar to those reported in double-diffusive convection flows. Our analysis demonstrates that this is, in fact, a transient regime. The third regime emerges at very large ${\varLambda }$, where temperature and concentration in the liquid undergo fully diffusive transport, leading to the elimination of convection flow due to the developing stably stratified layer. Notably, the diffusion regime can further be categorized into ‘melting’ or ‘dissolution’ patterns, each dominated by temperature or solute concentration, exhibiting distinct characteristics.

Furthermore, we conduct a quantitative analysis to elucidate the flow characteristics of the system. In the convection regime, both the thermal and solutal Nusselt numbers exhibit power-law scaling, approximately proportional to ${\sim }{Ra}_e^{1/4}$. The prefactors of these Nusselt numbers are thoroughly understood through the establishment of an energy and mass budget. As the system transitions into the diffusion regime, ${Nu}_T$ converges to $1$, while ${Nu}_C$ converges to $2.49$ for the extreme ‘melting’ situation, and $1$ for the extreme ‘dissolution’ scenario. Next, we delve into the interfacial values of temperature and concentration, revealing their dependence on specific flow regimes.

We then investigate the temporal development of the interface height $\bar {h}(t)$ across different regimes. In the convection and the ‘melting’-dominated diffusion regimes, we use the energy budget and the correlation ${Nu}_T(t) = \gamma _1\,{Ra}_e^{\gamma _2}(t)$ to derive that $\bar {h}(t) \sim t^{4/5}$ and $\bar {h}(t) \sim t^{1/2}$, respectively. Accordingly, the melting rate of the solid follows scaling laws $t^{-1/5}$ and $t^{-1/2}$ in these two regimes. In the ‘dissolution’-dominated diffusion regime, the interaction between the interfacial temperature and concentration necessitates solving a nonlinear ordinary differential equation to determine $\bar {h}(t)$. Additionally, in the convection regime, the roughness of the interface is observed to scale as $h_{{rms}}(t) \sim {Ra}_e^{1/3}$, consistent with findings reported for the one-component melting problem (Yang et al. 2023b). These results are crucial for developing models to predict $\bar {h}(t)$ in real-world scenarios. For example, we predict that the melting rate in the ‘melting’-dominated diffusion regime is $t^{-1/2}$, which aligns well with the ice-shelf melting results reported by Rosevear et al. (2022).

Despite the richness of the flow structures revealed in this study, there are avenues for more systematic investigations in the future. For example, exploring scenarios with a lighter solute, leading to a two-scalar RBC with a destabilizing solutal effect (Radko 2013), could provide further insights. Geophysically, considering scenarios such as the melt in a magma ocean where ${Pr}\rightarrow \infty$ while the liquid metal has ${Pr} \sim O(10^{-2})$, the impact of ${Pr}$ could lead to significant changes in flow structures (van der Poel et al. 2013). Additionally, despite the computational resources required, three-dimensional simulations would still be valuable in unveiling the three-dimensional structures of similar problems.