Theoretical and numerical studies on dissolution of horizontal salt body and pattern formation incorporated with a dynamic moving interface

Abstract

The dissolution of a horizontal salt body is theoretically and numerically investigated by considering the coupled dynamics of the dissolution reaction and the diffusive and convective mass transfer at the salt surface. To take into account the nature of the recessing salt surface and the dissolution-driven convection, a dynamic moving boundary condition coupled with a dissolution reaction are employed By employing a newly derived moving interface condition and parameters, a theoretical analysis predicts the onset of gravitational instability, suggesting scaling relationships between parameters to describe the onset of instability (in transport-dominant (Damkholer number $Da \to \infty $ and Rayleigh number $Ra > {10^5}$), onset time ${\tau _d}\sim R{a^{ - 2/3}}$ and ${\tau _d}\sim (RaDa)^{1/4}$ for reaction-dominant regimes $(Da \to 0)$). The scaling relationships on dissolution (the average height change, $\Delta {h_{avg}}\sim \sqrt \tau $ for the diffusion-dominant regime and $\Delta {h_{avg}}\sim R{a^{0.79/3}}{\tau ^{0.86}}$ for the convection-dominant regime) are also suggested to explain the pattern formation. Under a convective mass transfer condition, the concentration gradient of the solute developed during the recession of the salt surface results in an adverse density gradient, which causes gravitational instability motion and finally determines the formation of a dissolution pattern. In addition, two-dimensional and three-dimensional numerical simulations visualize dissolution-driven convective flow fields, the recession of the liquid-solid interface and pattern formation on the dissolving interface. Theoretical and numerical analyses of the dissolution process suggest that the control of gravitational instability is important to prevent or enhance dissolution pattern formation. This systematic parameter study provides a deep understanding of the effect of gravitational instability on convection-driven dissolution in a horizontal geometry.

1. Introduction

The dissolution of a solid in a liquid is ubiquitous in nature and industry such as erosion, rock fractures, melting, etc. The dissolution rate and periodic pattern formation of a soluble solid body are determined by the complex mass transfer and reaction rate (Wagner 1949; Cohen et al. 2016; Guérin et al. 2020; Hu et al. 2021; Wang et al. 2022). Specifically, when a salt is dissolving in a liquid, dynamic mass transport phenomena around the dissolving solid exist, changing the shape of the solid body of the salt. Huang et al. (2015) remarkably demonstrated the recession of a dissolving solid in a high-speed laminar flow and the effect of the solid body geometry on the dissolution process. Accordingly, to predict the dissolution rate and pattern formation successfully, several factors, including the flow and body geometry, must be considered in addition to the mass transfer on the dissolving surface (Kimura, Tanaka & Sukegawa 1990; Sullivan, Liu & Ecke 1996; Huang et al. 2015; Nakouzi, Goldstein & Steinbock 2015; Philippi et al. 2019; Cohen et al. 2020). The effects of gravitational instability on the dissolution rate and pattern formation under a stationary solvent condition have been studied, but only in a limited manner.

Nakouzi et al. (2015) showed that a vertical water-soluble cylinder object evolves into a paraboloid shape due to the inevitable density-driven flow in the vertical geometry, which experiences non-uniform flux. Later, Pegler & Davies Wykes (2020) analysed theoretically and experimentally the effect of stable laminar natural convective flow on the shape evolution of vertical conic solid bodies. In addition, Davies Wykes et al. (2018) experimentally examined the shape dynamics of a dissolving body in the presence of the buoyancy driven self-generated flows. In their experiments, for various quasi-two-dimensional and three-dimensional (3-D) axisymmetric solids, the dissolution of the underside is strongly affected by the gravitational instabilities generated by the unstable density profile near the surface and, therefore, on the underside of a partially dissolved wedge, small-scale roughness is observed. In the context of ultra-sharp pinnacle formation, Huang et al. (2020) experimentally and theoretically studied the effects of thermodynamics (Joule–Thompson effect) and hydrodynamics (natural convective flow) on the formation of pinnacles. According to their analysis, the hydrodynamic factor is important in the formation of ultra-sharp pinnacles.

In an inclined geometry, a parallel groove-type dissolution pattern forms at the surface of an initially flat soluble material due to a forced convective water flow (Guérin et al. 2020). Cohen et al. (2020) studied the effect of a natural convective flow on the dissolution rate of an inclined solid block and pattern formation on the dissolving surface, finding that parallel stripes form beyond a certain length from the block end. Over time, these stripes cross and evolve into scallops that propagate upstream. Meanwhile, in a horizontal geometry, the formation of a dissolution pattern is affected by gravitational instability in a different manner. Sullivan et al. (1996) experimentally visualized the dissolution of horizontal crystal bodies of NaCl, KBr and KCl in aqueous solutions. The dissolution rate and surface profile were compared with theoretical predictions in terms of Howard's (1966) boundary-layer instability model and Foster's (1968) linear stability. In an analysis with the boundary-layer instability model, the turbulent mass transfer rate was assumed to be closely related to the stability of the fluid layer. It was found that the wavelength of the dissolution pattern is significantly affected by the concentration gradient near the surface. Philippi et al. (2019) specifically identified the effects of dissolution-driven flow instability by defining the sequence of the dynamic transport phenomena around a solid body with three regimes: (i) one that was diffusive, (ii) one with the growth of instability and (iii) one that featured the emission of a plume, termed a quasi-stationary regime. Similar to the work by Sullivan et al. (1996), Philippi et al. (2019) introduced a boundary-layer instability model in the quasi-stationary regime and linked the stability of the liquid layer to the mass transfer rate. Under the assumption of high mass transfer rates in the turbulent natural convection regime (Sullivan et al. 1996) and the quasi-stationary regime (Philippi et al. 2019), they discussed the effects of the dissolution reaction on the surface on the onset of gravitational instability during dissolution and proposed several scaling relationships and compared these outcomes with experimental results. In addition, under the static assumption that the movement of the solid interface is much slower than the dissolution-driven convective motion, Philippi et al. (2019) and Berhanu et al. (2021) conducted nonlinear numerical simulations and a linear stability analysis of dissolution, respectively. Due to such a simplification, they could not simulate the temporal evolution of the interface position, i.e. pattern formation by dissolution-enhanced convective motion. Hence, none of the analyses considering realistic moving boundary conditions are applicable to simulate the dissolution process. Thus, while fluid flows are known to promote the dissolution of materials, such processes are poorly understood due to the coupled dynamics of the flow and the receding surface.

In this study, considering a moving interfacial boundary on a dissolving solid, the effects of dissolution on the onset of buoyancy-driven gravitational instability are coupled to simulate the dissolution of a horizontal salt body using linear stability theory and a numerical simulation. According to the dissolution regimes (i) and (ii) defined by Philippi et al. (2019), the effect of gravitational instability on the interfacial pattern formation on a dissolving surface is systemically studied. To simulate the growth of instability, a two-dimensional (2-D) nonlinear analysis is conducted using the commercial finite element method solver COMSOL Multiphysics (2019). As a case study, the dissolution of NaCl salt in water is simulated. Furthermore, the dissolution-driven interfacial pattern formation on a dissolving salt surface is visualized by fully 3-D numerical simulations. The effects of gravitational instability on the dissolution rate and pattern formation are demonstrated successfully. The present study will provide basic tools with which to predict the dissolution rate and to understand the dissolution-driven pattern formation process during the dissolution of a horizontal salt in a stationary solvent.

2. Governing equations and base fields

The system considered here is that of a horizontal salt body floating on top of a liquid with the initial liquid thickness of d and the initial concentration of solute in the liquid of ${C_r}$. The schematic diagram of the system is shown in figure 1. Initially, the liquid is bounded by a flat solid salt. As time goes on $(t \ge 0)$, the solid salt dissolves in liquid at the solid salt–liquid interface until ${C_r}$ is lower than the saturation concentration, ${C_{sat}}$. In the meantime, the position of the interface moves upward until ${C_r} = {C_{sat}}$. In a binary mixture such as aqueous solution of salt, a linear relation between the density $\rho $ and mass concentration of solute C in the solution has been assumed (Philippi et al. 2019). Under this assumption, the density of solution can be expressed as

(2.1)\begin{equation}\rho = {\rho _r} + \beta C,\end{equation}

where ${\rho _r}$ is the density of pure solvent and $\beta \{ = ({\rho _{sat}} - {\rho _r})/{C_{sat}}\} $ is the densification coefficient. The concentration gradient driven by the salt dissolution makes the system unstable and convective motion will begin in a certain time period if the dissolution rate of salt is faster than the diffusive transport rate of dissolved salt. Considering sufficiently dilute concentrations of solute, we use the Boussinesq approximation, which consists of treating the liquid density as constant in all terms in the equations of motions except the one in the gravity force (Chandrasekhar 1961). The flow is then incompressible. By assuming that the liquid is a Newtonian fluid and the variation of temperature of the liquid can be neglected (Sullivan et al. 1996), the governing equations of flow and concentration fields are expressed as

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

(2.3)\begin{gather}{\rho _r}\left\{ {\frac{\partial }{{\partial t}} + \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla }} \right\}\boldsymbol{U} =- \boldsymbol{\nabla }P + \boldsymbol{\nabla }\boldsymbol{\cdot }\{ \mu \boldsymbol{\nabla }\boldsymbol{U} + \mu {(\boldsymbol{\nabla }\boldsymbol{U})^T}\} + \beta C\boldsymbol{g},\end{gather}

(2.4a)\begin{gather}\frac{{\partial C}}{{\partial t}} =- \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{N},\end{gather}

(2.4b)\begin{gather}\boldsymbol{N} =- \mathrm{{\mathcal{D}}}\boldsymbol{\nabla }C + C\boldsymbol{U},\end{gather}

(2.5)\begin{gather}\mu = \mu _{r} \exp(\gamma (C - C_{r})).\end{gather}

Here, $\boldsymbol{U}( = \boldsymbol{i}U + \boldsymbol{j}V + \boldsymbol{k}W)$, P, $\boldsymbol{g}$, $\mu $, $C$, $\gamma$ and $\mathrm{{\mathcal{D}}}$ represent the velocity vector, the pressure, the gravitational acceleration, the viscosity, the mass concentration, the viscosity-variation parameter and the diffusivity, respectively. The variation of viscosity depending on concentration is modelled by (2.5), which has been widely used in related fields (Nijjer, Hewitt & Neufeld 2018; Kim et al. 2021). The given initial and boundary conditions are

(2.6a)\begin{gather}\boldsymbol{U} = 0\quad \textrm{and}\quad C = {C_r}\;\textrm{at}\;t = 0,\end{gather}

(2.6b)\begin{gather}\boldsymbol{U} = {\boldsymbol{U}_L}\quad \textrm{and}\quad C = {C_i},\;\textrm{at}\;Z = H(t,X,Y),\end{gather}

(2.6c)\begin{gather}\boldsymbol{U} = 0\quad \textrm{and}\quad \frac{{\partial C}}{{\partial Z}} = 0\;\textrm{at}\;Z = d,\end{gather}

where $H(t,X,Y)$ is the position of dissolution front, ${C_i}$ is the concentration at the solid salt–liquid interface and ${\boldsymbol{U}_L}$ is the velocity of the liquid phase at the interface, which will be discussed below.

Figure 1. Schematic diagram of the system studied here. Initially, the interface between the solid salt and the solution is located at $Z = 0$. As the salt dissolves into the solution, the interface position $H(t)$ moves upward and the solute (dissolved salt) induces buoyancy-driven instability.

At the salt–liquid interface, the following kinematic condition can be derived from mass conservation (Leal 2007):

(2.7)\begin{equation}{\rho _S}({\boldsymbol{U}_s} - {\boldsymbol{U}_i})\boldsymbol{\cdot }\boldsymbol{n} = {\rho _L}({\boldsymbol{U}_L} - {\boldsymbol{U}_i})\boldsymbol{\cdot }\boldsymbol{n}\;\textrm{at}\;Z = H(t,X,Y),\end{equation}

where the subscript ‘S’ represents the properties of the solid salt, ${\boldsymbol{U}_i}$ is the interface moving velocity, $(\partial H/\partial t)$ is the normal dissolution rate and $\boldsymbol{n}$ is the unit normal vector. Since there is no motion within the solid salt phase, i.e. ${\boldsymbol{U}_s} = 0$, we can express ${\boldsymbol{U}_L}$ as

(2.8)\begin{equation}{\boldsymbol{U}_L} =- \frac{{\partial H}}{{\partial t}}(\chi - 1)\boldsymbol{n}\;\textrm{at}\;Z = H(t,X,Y),\end{equation}

where $\chi $ is the expansion factor, which is the ratio of the specific volume of the solute occupied in the liquid phase and that of the solid salt phase (Philippi et al. 2019; Cohen et al. 2020). At the moment, the density of liquid at the interface cannot be defined, however, it will be deduced from known quantities. Because the density of liquid at the interface has the same meaning as the interfacial concentration, from (2.4), ${\rho _L} = {\rho _r} + \beta {C_i} = {C_i}$. According to this relation, ${\rho _L}( = {C_i}) = {\rho _r}/(1 - \beta )$ and therefore, $\chi ( = {\rho _S}/{\rho _L}) = {\rho _S}/{\rho _r}(1 - \beta )$. The dissolution proceeds through the reaction at the surface and the mass transfer from the dissolving surface to bulk solution. Then, the dissolution rate can be expressed as (Philippi et al. 2019; Cohen et al. 2020)

(2.9a)\begin{gather}- {\rho _S}\frac{{\partial H}}{{\partial t}} = {r_{rxn}} = \alpha ({C_{sat}} - {C_i})\;\textrm{at}\;Z = H(t,X,Y),\end{gather}

(2.9b)\begin{gather}- {\rho _S}\frac{{\partial H}}{{\partial t}} = {r_{mt}} = \left( { - \frac{{\partial H}}{{\partial t}}{C_i}\boldsymbol{n} - \mathrm{{\mathcal{D}}}\boldsymbol{\nabla }C + \boldsymbol{U}{C_i}} \right)\boldsymbol{\cdot }\boldsymbol{n}\;\textrm{at}\;Z = H(t,X,Y),\end{gather}

where ${r_{rxn}}$ and ${r_{mt}}$ are the dissolution reaction rate and the mass transfer rate, respectively. The first term of the right-hand side of (2.9b) represents the convective transport due to the movement of the interface and the second term of right-hand side of (2.9b) is the diffusive and the convective transport in the liquid phase. Therefore, for a mass transfer-controlled dissolution case, such as NaCl dissolution in aqueous solution (Wagner 1949), the dissolution front can be determined by the following material balance (Cohen et al. 2020):

(2.10)\begin{equation}- {\rho _S}\frac{{\partial H}}{{\partial t}} = \left\{ { - \mathrm{{\mathcal{D}}}\boldsymbol{\nabla }C + \left( { - \frac{{\partial H}}{{\partial t}}} \right)\chi {C_{sat}}\boldsymbol{n}} \right\}\boldsymbol{\cdot }\boldsymbol{n}\;\textrm{at}\;Z = H(t,X,Y).\end{equation}

The important parameters governing the present system are the Schmidt number Sc, the Rayleigh number Ra and the Damkhöler number Da, which are defined as

(2.11a–c)\begin{equation}Sc = \frac{\upsilon }{D},\quad Ra = \frac{{g\beta {C_{sat}}{d^3}}}{{\mathrm{{\mathcal{D}}}\mu }}\quad \textrm{and}\quad Da = \frac{{\alpha d}}{D},\end{equation}

where $\upsilon $ is the kinematic viscosity of solvent. The Damkhöler number represents the ratio between the reaction rate and the diffusive mass transfer rate.

Employing d, ${d^2}/\mathrm{{\mathcal{D}}}$ and $({C_{sat}} - {C_r})$ as length, time and concentration scaling factors, the above governing equations (2.1)–(2.5) are non-dimensionalized as

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

(2.13)\begin{gather}\frac{1}{{Sc}}\left\{ {\frac{\partial }{{\partial \tau }} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }} \right\}\boldsymbol{u} =- \boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }\{ \bar{\mu }\boldsymbol{\nabla }\boldsymbol{u} + \bar{\mu }{(\boldsymbol{\nabla }\boldsymbol{u})^T}\} + Ra\,c\boldsymbol{k},\end{gather}

(2.14)\begin{gather}\left( {\frac{\partial }{{\partial \tau }} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }} \right)c = {\nabla ^2}c,\end{gather}

where $\bar{\mu }( = \mu /{\mu _r}) = \exp (\varGamma c)$ and $\varGamma = \gamma ({C_{sat}} - {C_r})$. The initial and boundary conditions (2.6) become

(2.15a)\begin{gather}\boldsymbol{u} = 0\quad \textrm{and}\quad c = 0\;\textrm{at}\;\tau = 0,\end{gather}

(2.15b)\begin{gather}\boldsymbol{u} = (1 - \chi )\left( {\frac{{\partial h}}{{\partial \tau }}} \right)\boldsymbol{n}\quad \textrm{and}\quad c = \frac{{{C_i}}}{{{C_{sat}} - {C_r}}}\;\textrm{at}\;z = h(\tau ,x,y),\end{gather}

(2.15c)\begin{gather}\boldsymbol{u} = 0\quad \textrm{and}\quad \frac{{\partial c}}{{\partial z}} = 0\;\textrm{at}\;z = 1.\end{gather}

From (2.6b) and (2.8), the solid–liquid interface position, $h(\tau )$, can be determined from (2.9) as

(2.16a,b)\begin{equation}\frac{{\partial h}}{{\partial \tau }} =- Da\frac{{\Delta C}}{{{\rho _S}}}\left( {\frac{{{C_{sat}}}}{{\Delta C}} - c} \right)\quad \textrm{and}\quad \frac{{\partial h}}{{\partial \tau }} = \frac{{\Delta C/{\rho _S}}}{{1 - \chi (\Delta C/{\rho _S})c}}\boldsymbol{\nabla }c\boldsymbol{\cdot }\boldsymbol{n}\;\textrm{at}\;z = h(\tau ,x,y).\end{equation}

Then, the following auxiliary condition can be derived:

(2.17)\begin{equation}\boldsymbol{\nabla }c\boldsymbol{\cdot }\boldsymbol{n} =- Da\textrm{(}1 - c\textrm{)\{ }1 - \chi (\mathrm{\Delta }C/{\rho _s})c\textrm{\} }.\end{equation}

For mathematical simplicity, from now on, we assume ${C_r} = 0$, i.e. initially the liquid is free from the solute, then in (2.16a) and (2.16b) ${C_{sat}}/\Delta C = 1$ and $\Delta C = {C_{sat}}$. In the mass transfer systems, $Sc \gg 1$ has been assumed as usual.

Before the onset of instability, the base-concentration profile can be governed by

(2.18a)\begin{gather}\frac{{\partial {c_0}}}{{\partial \tau }} + {w_0}\frac{{\partial {c_0}}}{{\partial z}} = \frac{{{\partial ^2}{c_0}}}{{\partial {z^2}}},\end{gather}

(2.18b)\begin{gather}{w_0} = (1 - \chi )\left( {\frac{{\partial {h_0}}}{{\partial \tau }}} \right),\end{gather}

under the following initial and boundary conditions:

(2.19a)\begin{gather}{c_0}(\tau ,z) = 0\;\textrm{at}\;\tau = 0,\end{gather}

(2.19b)\begin{gather}{c_0} = 1 + \frac{1}{{Da}}\frac{1}{{1 - \chi ({C_{sat}}/{\rho _s}){c_0}}}\frac{{\partial {c_0}}}{{\partial z}}\;\textrm{at}\;z = {h_0}(\tau ),\end{gather}

(2.19c)\begin{gather}\frac{{\partial {c_0}}}{{\partial z}} = 0\;\textrm{at}\;z = 1,\end{gather}

(2.19d)\begin{gather}\frac{{\partial {h_0}}}{{\partial \tau }} = \frac{{{C_{sat}}/{\rho _s}}}{{1 - \chi ({C_{sat}}/{\rho _s}){c_0}}}{\left. {\frac{{\partial {c_0}}}{{\partial z}}} \right|_{{h_0}}}.\end{gather}

Under the static assumption, i.e. ${h_0} = 0$ and ${w_0} = 0$, Philippi et al. (2019) and Berhanu et al. (2021) suggested the following base-concentration field:

(2.20)\begin{equation}{c_0} = erfc\left( {\frac{z}{{2\sqrt \tau }}} \right) - \exp (Daz + D{a^2}\tau )erfc\left( {\frac{z}{{2\sqrt \tau }} + Da\sqrt \tau } \right),\end{equation}

under the boundary condition

(2.21)\begin{equation}\boldsymbol{\nabla }c\boldsymbol{\cdot }\boldsymbol{n} =- Da(1 - c)\;\textrm{at}\;z = 0.\end{equation}

It should be kept in mind that, for the extreme cases of $\chi \to 1$ and ${C_{sat}}/{\rho _s} \to 0$, Philippi et al.'s (2019) and Berhanu et al.'s (2021) static boundary condition (2.21) can be deduced form the present boundary conditions (2.16) and (2.17). Furthermore, for the limiting case of $D{a^\ast }( = Da\sqrt \tau ) \to \infty $, the above solution (2.20) is reduced to

(2.22)\begin{equation}{c_0} = erfc\left( {\frac{z}{{2\sqrt \tau }}} \right).\end{equation}

Because the present base-concentration field, (2.18) and (2.19), is a complex function of $D{a^\ast }$, $({C_{sat}}/{\rho _S})$ and $\chi $, we will consider two limiting cases of $Da \to \infty $ (transport-controlled system) and $Da \to 0$ (reaction-controlled system) first in the next sections.

2.1. For the limiting case of $Da \to \infty $

For this transport-controlled case, from (2.16a) $c \to 1$ at $z = h(\tau ,x,y)$ is obtained, i.e. ${C_i} \to {C_{sat}}$. In this case, the boundary condition (2.15b) and the interface movement equation (2.16) can be rewritten as

(2.23a,b)\begin{equation}\frac{{\partial h}}{{\partial \tau }} = \frac{1}{{1 - \chi ({C_{sat}}/{\rho _S})}}\boldsymbol{\nabla }c\boldsymbol{\cdot }\boldsymbol{n}\quad \textrm{and}\quad c = 1,\;\textrm{at}\;z = h(\tau ,x,y).\end{equation}

Since the above (2.18) and (2.19) are quite similar to the Stefan problem (Carslaw & Jaeger 1959), under the assumption of a deep pool, i.e. $\sqrt \tau \ll 1$, the above concentration field of the liquid phase can be represented as (see § 11.2 of Carslaw & Jaeger (1959))

(2.24a–c)\begin{align}{c_0} &= Aerfc\left( {\frac{z}{{2\sqrt \tau }} - {w_0}\sqrt \tau } \right),\quad {h_0}(\tau ) =- \eta \sqrt \tau \quad \textrm{and}\notag\\ {w_0} &= (1 - \chi )\left( {\frac{{\partial {h_0}}}{{\partial \tau }}} \right) =- \frac{{(1 - \chi )\eta }}{{2\sqrt \tau }},\end{align}

where $erfc(x) = 1 - erf(x)$ and $erf(x)$ is the error function. From (2.16a) and (2.17), the undetermined constants A and $\eta $ can be obtained by solving the following nonlinear simultaneous equations:

(2.25a) \begin{gather}Aerfc\left( { - \frac{{\chi \eta }}{2}} \right) = 1 - \frac{1}{{D{a^\ast }}}\frac{A}{{1 - \chi ({C_{sat}}/{\rho _s})}}\frac{1}{{\sqrt {\rm \pi} }}\exp \left\{ { - {{\left( { - \frac{{\chi \eta }}{2}} \right)}^2}} \right\},\end{gather}

(2.25b)\begin{gather}\frac{\eta }{2}\sqrt {\rm \pi} \exp \left\{ {{{\left( { - \frac{{\chi \eta }}{2}} \right)}^2}} \right\} = \frac{{A({C_{sat}}/{\rho _s})}}{{1 - \chi ({C_{sat}}/{\rho _s})Aerfc( - \chi \eta /2)}},\end{gather}

where $D{a^\ast } = Da\sqrt \tau $.

For the transport-controlled case, from (2.21), the base-concentration field can be given as

(2.26a)\begin{gather}{c_0} = \frac{1}{{erfc( - \chi \eta /2)}}erfc\left\{ {\frac{z}{{2\sqrt \tau }} + \frac{{\eta (1 - \chi )}}{2}} \right\},\end{gather}

(2.26b)\begin{gather}\frac{{({C_{sat}}/{\rho _S})}}{{1 - \chi ({C_{sat}}/{\rho _S})}} = \frac{\eta }{2}\sqrt {\rm \pi} \exp \left\{ {{{\left( {\frac{{ - \chi \eta }}{2}} \right)}^2}} \right\}erfc\left( { - \frac{{\chi \eta }}{2}} \right),\end{gather}

except for the singular limit of $\tau = 0$.

For the limiting case of $\chi = 1$, where the volume expansion can be neglected, Verhaeghe et al. (2005) suggested the following relation:

(2.27a)\begin{gather}{c_0} = \frac{1}{{erfc( - \eta /2)}}erfc\left\{ {\frac{z}{{2\sqrt \tau }}} \right\},\end{gather}

(2.27b)\begin{gather}\frac{{({C_{sat}}/{\rho _S})}}{{1 - ({C_{sat}}/{\rho _S})}} = \frac{\eta }{2}\sqrt {\rm \pi} \exp \left( {\frac{{{\eta^2}}}{4}} \right)erfc\left( { - \frac{\eta }{2}} \right).\end{gather}

For another limiting of $\chi \to 0$, where ${\rho _L} \gg {\rho _S}$, (2.26) can be reduced as

(2.28a,b)\begin{equation}{c_0} = erfc\left( {\frac{z}{{2\sqrt \tau }} + \frac{\eta }{2}} \right)\quad \textrm{and}\quad \eta = \frac{2}{{\sqrt {\rm \pi} }}\left( {\frac{{{C_{sat}}}}{{{\rho_S}}}} \right).\end{equation}

As expected, the concentration field depends strongly on $({C_{sat}}/{\rho _S})$ and $\chi $. In figure 2(a), $\eta $ is plotted as functions of $({C_{sat}}/{\rho _S})$ and $\chi $. For the extreme case of $({C_{sat}}/{\rho _S}) \to 0$, the above concentration fields (2.26) are reduced to (2.22), since $\eta \to 0$ as $({C_{sat}}/{\rho _S}) \to 0$. For the limiting case of $Da \to \infty $ and $({C_{sat}}/{\rho _S}) \to 0$, the present base field is the same as (2.20). This means that Berhanu et al.'s (2021) analysis has its own limit, because their base-concentration profile, (2.20), is quite different from the present ones, (2.26). In figures 2(b) and 2(c), the present solutions are compared with the previous one, (2.24), for the limiting cases of $D{a^\ast } \to \infty $. Unlike the present study, the previous solution (2.24) is independent of $\chi $. Except for the limiting case of $({C_{sat}}/{\rho _S}) \to 0$, the present base-concentration fields are quite different from the previous solution under the static assumption.

Figure 2. Comparison of base-concentration fields of various systems. (a) Prediction of interface movement parameter $\eta $ of transport-controlled systems by solving (2.25a) and (2.25b) and comparison of base-concentration fields for various systems: (b) case of $D{a^\ast } \to \infty $ and $\chi = 1$ (solution of (2.27)); (c) case of $D{a^\ast } \to \infty $ and $\chi = 0$ (solution of (2.28)); and (d) case of $D{a^\ast } \to 0$.

2.2. For the limiting case of $Da \to 0$

For another limiting case of $Da \to 0$, i.e. reaction-controlled system, the auxiliary boundary condition, (2.17), can be rewritten as $\boldsymbol{\nabla }(c/Da)\boldsymbol{\cdot }\boldsymbol{n} =- \{ 1 - Da(c/Da)\} \{{1 - Da\chi (\Delta C/{\rho_s})(c/Da)} \}$. Since, from (2.16a) $\partial h/\partial \tau \to 0$ as $Da \to 0$ can be obtained, $({C_i}/{C_{sat}}) \to 0$ is assumed for this very slow dissolution rate system. In this case, the governing equations (2.13) and (2.14), and the boundary condition (2.15b) can be rewritten as

(2.29)\begin{gather}\frac{1}{{Sc}}\left\{ {\frac{\partial }{{\partial \tau }} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }} \right\}\boldsymbol{u} =- \boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }\{ \bar{\mu }\boldsymbol{\nabla }\boldsymbol{u} + \bar{\mu }{(\boldsymbol{\nabla }\boldsymbol{u})^T}\} + R{a_D}\tilde{c}\boldsymbol{k},\end{gather}

(2.30)\begin{gather}\left( {\frac{\partial }{{\partial \tau }} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }} \right)\tilde{c} = {\nabla ^2}\tilde{c},\end{gather}

under the following boundary conditions:

(2.31)\begin{equation}\boldsymbol{u} = 0\quad \textrm{and}\quad \frac{{\partial \tilde{c}}}{{\partial z}} =- 1\;\textrm{at}\;z = 0,\end{equation}

where $\tilde{c} = c/Da$ and $R{a_D}( = Ra \times Da) = g\beta \alpha {C_{sat}}{d^4}/{D^2}\upsilon $.

In this limiting case of $Da \to 0$, the base-concentration profile can be obtained by solving the following diffusion equation:

(2.32)\begin{equation}\frac{{\partial {{\tilde{c}}_0}}}{{\partial \tau }} = \frac{{{\partial ^2}{{\tilde{c}}_0}}}{{\partial {z^2}}},\end{equation}

under the initial and boundary conditions

(2.33a–c)\begin{equation}{\tilde{c}_0}(0,z) = 0,\quad \frac{{\partial {{\tilde{c}}_0}}}{{\partial z}}(\tau ,0) =- 1\quad \textrm{and}\quad \frac{{\partial {{\tilde{c}}_0}}}{{\partial z}}(\tau ,1) = 0,\end{equation}

under the assumption of a deep pool, i.e. $\sqrt \tau \ll 1$. By using the Laplace transform method, (2.32) and (2.33) can be solved as

(2.34)\begin{equation}{c_0}/D{a^\ast } = {\tilde{c}_0}(\tau ,z)/\sqrt \tau = 2ierfc\left( {\frac{z}{{2\sqrt \tau }}} \right),\end{equation}

where $ierfc(x) = 1/\sqrt {\rm \pi} \exp ( - {x^2}) - xerfc(x)$. The above solution is featured in figure 2(d). For this limiting case, Philippi et al.'s (2019) and Berhanu et al.'s (2021) base-concentration field, (2.20), yields the following unexpected result:

(2.35)\begin{equation}{c_0} = 0.\end{equation}

Then, for limiting case of $Da \to 0$, further analysis is not possible using Philippi et al.'s (2019) and Berhanu et al.'s (2021) base-concentration field.

3. Linear stability analysis

3.1. Stability equations

Under linear stability theory, infinitesimal disturbances caused by incipient convective motion at the dimensionless critical time ${\tau _c}$ can be formulated in terms of the concentration component ${c_1}$ and the vertical velocity component ${w_1}$ by linearizing equations (2.12)–(2.14), and then taking double curl on the linearized equations of (2.12) and (2.13)

(3.1)\begin{align}\frac{1}{{Sc}}\left( {\frac{\partial }{{\partial \tau }} + {w_0}\frac{\partial }{{\partial z}}} \right){\nabla ^2}{w_1} &= {\bar{\mu }_0}{\nabla ^4}{w_1} + 2\frac{{\partial {{\bar{\mu }}_0}}}{{\partial z}}{\nabla ^2}\left( {\frac{{\partial {w_1}}}{{\partial z}}} \right) + \frac{{{\partial ^2}{{\bar{\mu }}_0}}}{{\partial {z^2}}}\left( {{\nabla^2}{w_1} - 2\frac{{{\partial^2}{w_1}}}{{\partial {z^2}}}} \right)\notag\\ &\quad + Ra\nabla _1^2{c_1},\end{align}

(3.2)\begin{gather}\frac{{\partial {c_1}}}{{\partial \tau }} + {w_0}\frac{{\partial {c_1}}}{{\partial z}} + {w_1}\frac{{\partial {c_0}}}{{\partial z}} = {\nabla ^2}{c_1},\end{gather}

where ${\nabla ^2} = {\partial ^2}/\partial {z^2} + \nabla _1^2$ and $\nabla _1^2 = {\partial ^2}/\partial {x^2} + {\partial ^2}/\partial {y^2}$. The proper boundary conditions are

(3.3a)\begin{gather}{w_1} = (1 - \chi )\left( {\frac{{\partial {h_1}}}{{\partial \tau }}} \right)\quad \textrm{and}\quad \frac{{\partial {h_1}}}{{\partial \tau }} = Da\frac{{{C_{sat}}}}{{{\rho _S}}}{c_1}\;\textrm{at}\;z = {h_0}(\tau ),\end{gather}

(3.3b)\begin{gather}\frac{{\partial {c_1}}}{{\partial z}} = Da{c_1}\{ 1 - \chi ({C_{sat}}/{\rho _S}){c_0}\} + Da(1 - {c_0})\chi ({C_{sat}}/{\rho _S}){c_1}\;\textrm{at}\;z = {h_0}(\tau ),\end{gather}

(3.3c)\begin{gather}{w_1} = \frac{{\partial {w_1}}}{{\partial z}} = \frac{{\partial {c_1}}}{{\partial z}} = 0\;\textrm{at}\;z = 1.\end{gather}

The disturbances are assumed to exhibit the horizontal periodicity, then the following Fourier mode is employed:

(3.4)\begin{equation}[{w_1}(\tau ,x,y,z),{c_1}(\tau ,x,y,z)] = [{w_1}(\tau ,z),{c_1}(\tau ,z)]\textrm{exp\{ }i({k_x}x + {k_y}y)\textrm{\} },\end{equation}

where ‘$i$’ is the imaginary number and ${k_x}$ and ${k_y}$ are the horizontal wavenumbers in the $x$- and $y$-directions. Using this Fourier mode analysis, the horizontal Laplacian becomes $\nabla _1^2 =- {k^2}$, where the horizontal wavenumber ‘$k$’ has the relation of $k = \sqrt {k_x^2 + k_y^2} $.

For the mass transfer system investigated in this study, the inertia terms have been neglected by assuming $Sc \to \infty $. Moreover, for the case of large $Ra$, the incipient convective motion is confined within the narrow penetration region, $\Delta \sim \sqrt {Dt} $ (Kim, Park & Choi 2002). In this case, the dominant operator of the global $(\tau ,z)$-domain – ${\partial ^2}/\partial {z^2}$ – does not have eigenfunctions that are localized around the base-concentration front (Riaz et al. 2006). Therefore, we transformed the disturbance equations such that the eigenfunctions associated with the streamwise diffusion operator are localized around the base-concentration front. Following a coordinate transformation to the similarity variable of the base state $\zeta ( = z/\sqrt \tau )$, for the case of $Sc \gg 1$ and $\sqrt \tau \ll 1$, the stability equations (3.1)–(3.4) can be reformulated using the relations of $\partial /\partial \tau {|_z} = \partial /\partial \tau {|_\zeta } + (\partial /\partial \zeta )(\partial \zeta /\partial \tau )$ and $\partial \zeta /\partial \tau =- \zeta /(2\tau )$ as

(3.5)\begin{gather}\left\{ {{{\bar{\mu }}_0}\left( {\frac{{{\partial^2}{w_1}}}{{\partial {\zeta^2}}} - {k^{{\ast} 2}}} \right) + 2\frac{{\partial {{\bar{\mu }}_0}}}{{\partial \zeta }}} \right\}\left( {\frac{{{\partial^2}{w_1}}}{{\partial {\zeta^2}}} - {k^{{\ast} 2}}} \right)w_1^\ast- \frac{{{\partial ^2}{{\bar{\mu }}_0}}}{{\partial {\zeta ^2}}}\left( {\frac{{{\partial^2}{w_1}}}{{\partial {\zeta^2}}} + {k^{{\ast} 2}}} \right)w_1^\ast= {k^{{\ast} 2}}{c_1},\end{gather}

(3.6)\begin{gather}\tau \frac{{\partial {c_1}}}{{\partial \tau }} + {w_0}\frac{{\partial {c_1}}}{{\partial \zeta }} + R{a^\ast }{w_1}\frac{{\partial {c_0}}}{{\partial \zeta }} = ({\mathrm{{\mathcal{L}}}_\zeta } - {k^{{\ast} 2}}){c_1},\end{gather}

under the following boundary conditions:

(3.7a)\begin{gather}w_1^\ast= \frac{{(1 - \chi )}}{{R{a^\ast }}}D{a^\ast }\left( {\frac{{{C_{sat}}}}{{{\rho_S}}}} \right){c_1}\quad \textrm{and}\quad \frac{{\partial w_1^\ast }}{{\partial \zeta }} = 0\;\textrm{at}\;\zeta = {h_0}/\sqrt \tau ,\end{gather}

(3.7b)\begin{gather}\frac{{\partial {c_1}}}{{\partial \zeta }} = D{a^\ast }{c_1}\{ 1 - \chi ({C_{sat}}/{\rho _S}){c_0}\} + D{a^\ast }(1 - {c_0})\chi ({C_{sat}}/{\rho _S}){c_1}\;\textrm{at}\;\zeta = {h_0}/\sqrt \tau ,\end{gather}

(3.7c)\begin{gather}w_1^\ast \to 0,\quad \; \frac{{\partial w_1^\ast }}{{\partial \zeta }} \to 0\quad \textrm{and}\quad {c_1} \to 0\;\textrm{as}\;\zeta \to \infty ,\end{gather}

where ${k^\ast } = k\sqrt \tau $, $w_1^\ast= {w_1}/(Ra\tau )$, $R{a^\ast } = Ra{\tau ^{3/2}}$ and $D{a^\ast } = Da{\tau ^{1/2}}$. It should be kept in mind that the self-similarity applies only to the base-concentration field, and the amplitude and the spatial structure of disturbances are time dependent. This kind coordinate transform has been used in related fields (Ben, Demekhin & Chang 2002; Riaz et al. 2006). The newly introduced linear operator ${\mathrm{{\mathcal{L}}}_\zeta }$ is defined as

(3.8)\begin{equation}{\mathrm{{\mathcal{L}}}_\zeta } = \frac{{{\partial ^2}}}{{\partial {\zeta ^2}}} + \frac{\zeta }{2}\frac{\partial }{{\partial \zeta }}.\end{equation}

It should be noted that the present stability equations (3.5)–(3.7) in the $(\tau ,\zeta )$-domain are mathematically identical to (3.1)–(3.3) in the global $(\tau ,z)$-domain except for the singular limit of $\tau = 0$. However, unlike the trigonometric functions, which are the eigenfunctions of the linear operator in the global $(\tau ,z)$-domain, the eigenfunctions of the present linear operator ${\mathrm{{\mathcal{L}}}_\zeta }$ are localized around the base-concentration front. Due to this fact, further analytic approaches are possible in the present $(\tau ,\zeta )$-domain, which will be discussed in the following section.

3.2. Solution method

As discussed previously, we obtained the base-concentration fields only for the limiting cases of $D{a^\ast } \to \infty$ and $D{a^\ast } \to 0$. Therefore, these limiting cases will be solved first.

3.2.1. For the limiting case of $D{a^\ast } \to \infty$

For the simplest system with no vertical velocity at the dissolving interface and no viscosity change with concentration variation, i.e. $\chi = 1$, $\eta = 0$ and ${\bar{\mu }_0} = 1$, a fully analytic approach is possible. This is one of the advantages of the similarity transform which was introduced in the previous section.

For the transport-controlled system, from (2.22), the concentration gradient can be rewritten as

(3.9)\begin{equation}\frac{{\textrm{d}{c_0}}}{{\textrm{d}\zeta }} =- \frac{1}{{erfc( - \chi \eta /2)}}\frac{1}{{\sqrt {\rm \pi} }}\exp \left[ { - {{\left\{ {\frac{\zeta }{2} + \frac{{\eta (1 - \chi )}}{2}} \right\}}^2}} \right].\end{equation}

Moreover, the boundary conditions (3.3a) and (3.3b) can be reduced as

(3.10)\begin{equation}w_1^\ast= \frac{{(1 - \chi )}}{{R{a^\ast }}}\frac{{({C_{sat}}/{\rho _S})}}{{1 - \chi ({C_{sat}}/{\rho _S})}}\frac{{\partial {c_1}}}{{\partial \zeta }}\quad \textrm{and}\quad \frac{{\partial w_1^\ast }}{{\partial \zeta }} = {c_1} = 0\;\textrm{at}\;\zeta =- \eta .\end{equation}

It should be noted that $w_1^\ast= 0$ at $\zeta =- \eta $ for the case of $\chi = 1$. Furthermore, for the limiting case of $\chi = 1$ and $({C_{sat}}/{\rho _S}) \to 0$, the present stability problem degenerates to the transient Rayleigh–Bénard problem (Kim et al. 2002), where

(3.11)\begin{gather}\frac{{\textrm{d}{c_0}}}{{\textrm{d}\zeta }} = \frac{1}{{\sqrt {\rm \pi} }}\exp \left( { - \frac{{{\zeta^2}}}{4}} \right)\quad \textrm{and}\quad {h_0}(\tau ) = 0,\end{gather}

(3.12)\begin{gather}w_1^\ast= \frac{{\partial w_1^\ast }}{{\partial \zeta }} = {c_1} = 0\;\textrm{at}\;\zeta = 0.\end{gather}

For the case of $\chi = 1$, $({C_{sat}}/{\rho _S}) \to 0$ and ${\bar{\mu }_0} = 1$, using the generalized Fourier series we can express ${c_1}$ as

(3.13)\begin{equation}{c_1}(\tau ,\zeta ) = \sum\limits_{i = 1}^\infty {{A_i}(\tau ){\kappa _i}{\phi _i}(\zeta )} ,\end{equation}

where ${\kappa _i}$ is the normalization factor, $\{ {\phi _i}(\zeta )\} $ is the set of eigenfunctions of the operator ${\mathrm{{\mathcal{L}}}_\zeta }$ and they should satisfy the following Sturm–Liouville boundary value problem:

(3.14)\begin{equation}{\mathrm{{\mathcal{L}}}_\zeta }{\phi _i} = \frac{{{\textrm{d}^2}{\phi _i}}}{{\textrm{d}{\zeta ^2}}} + \frac{\zeta }{2}\frac{{\textrm{d}{\phi _i}}}{{\textrm{d}\zeta }} =- {\lambda _i}{\phi _i},\end{equation}

under the boundary conditions

(3.15)\begin{equation}{\phi _i}(0) = 0\quad \textrm{and}\quad \frac{{\textrm{d}{\phi _i}}}{{\textrm{d}\zeta }} \to 0\;\textrm{as}\;\zeta \to \infty .\end{equation}

Then, from (3.5) and (3.13), $w_1^\ast $ is expressed as

(3.16)\begin{equation}w_1^\ast (\tau ,\zeta ) = \sum\limits_{i = 1}^\infty {{A_i}(\tau ){\kappa _i}{\psi _i}({k^\ast },\zeta )} ,\end{equation}

where ${\psi _i}$ values can be obtained by solving

(3.17)\begin{equation}\left( {\frac{{{\textrm{d}^2}}}{{\textrm{d}{\zeta^2}}} - {k^\ast }^2} \right)\left\{ {\left( {\frac{{{\textrm{d}^2}}}{{\textrm{d}{\zeta^2}}} - {k^\ast }^2} \right){\psi_i}} \right\} = {k^{{\ast} 2}}{\phi _i},\end{equation}

under the following boundary conditions:

(3.18a)\begin{gather}{\psi _i} = \frac{{\textrm{d}{\psi _i}}}{{\textrm{d}\zeta }} = 0\;\textrm{at}\;\zeta = 0,\end{gather}

(3.18b)\begin{gather}{\psi _i} \to 0\quad \textrm{and}\quad \frac{{\textrm{d}{\psi _i}}}{{\textrm{d}\zeta }} \to 0\;\textrm{as}\;\zeta \to \infty .\end{gather}

The analytic solutions of (3.13)–(3.18) are summarized in Appendix A.

Substituting $c_1^\ast $ and $w_1^\ast $ into (3.6) and following the orthogonalization procedure, the stability equations are reduced to the following matrix form:

(3.19a)\begin{equation}\tau \frac{{\textrm{d}\boldsymbol{a}}}{{\textrm{d}\tau }} = \boldsymbol{Ba},\end{equation}

where

(3.19b)\begin{gather}{B_{ij}} =- ({\lambda _i} + {k^{{\ast} 2}}){\delta _{ij}} + R{a^\ast }\frac{1}{{\sqrt {\rm \pi} }}{H_{ij}},\end{gather}

(3.19c)\begin{gather}{H_{ij}} = \sqrt {\rm \pi} \int_0^\infty {{\phi _i}(\zeta )} \frac{{\textrm{d}{c_0}}}{{\textrm{d}\zeta }}{\psi _j}(\zeta )\exp \left( {\frac{{{\zeta^2}}}{4}} \right)\,\textrm{d}\zeta = \int_0^\infty {{\phi _i}(\zeta ){\psi _j}(\zeta )\,\textrm{d}\zeta } ,\end{gather}

(3.19d)\begin{gather}\boldsymbol{a} = {[{A_1},{A_2},{A_3}, \ldots ]^\textrm{T}},\end{gather}

for $i,j = 1,2, \ldots $. In (3.19c), $\textrm{exp(}{\zeta ^2}/4\textrm{)}$ is the weighting function of the Sturm–Liouville equation (3.14) and the variation of base field, $\textrm{d}{c_0}/\textrm{d}\zeta $, is given by (3.11). Here, it is stressed that the partial differential equations (3.5)–(3.7) are reduced to a set of ordinary differential equation (3.17), without any spatial discretization. After the integration by parts, the following relation is obtained:

(3.20) \begin{align} {H_{ij}} & = \int_0^\infty {{\phi _i}{\psi _j}\,\textrm{d}\zeta } = \dfrac{1}{{{k^{{\ast} 2}}}}\int_0^\infty {\left( {\dfrac{{{\textrm{d}^2}}}{{\textrm{d}{\zeta^2}}} - {k^{{\ast} 2}}} \right)} \left\{ {\left( {\dfrac{{{\textrm{d}^2}}}{{\textrm{d}{\zeta^2}}} - {k^{{\ast} 2}}} \right){\psi_i}} \right\}{\psi _j}\,\textrm{d}\zeta \notag\\ & = \dfrac{1}{{{k^{{\ast} 2}}}}\int_0^\infty {\left\{ {\left( {\dfrac{{{\textrm{d}^2}}}{{\textrm{d}{\zeta^2}}} - {k^{{\ast} 2}}} \right){\psi_i}} \right\}} \left\{ {\left( {\dfrac{{{\textrm{d}^2}}}{{\textrm{d}{\zeta^2}}} - {k^{{\ast} 2}}} \right){\psi_j}} \right\}\textrm{d}\zeta = {H_{ji}}. \end{align}

This means that the characteristic matrix $\boldsymbol{B}$ is normal, i.e. $\boldsymbol{B} = {\boldsymbol{B}^{\boldsymbol{T}}}$.

To trace the growth of the disturbance, the norm of the disturbance ${\|c_1\|}$ is defined as

(3.21)\begin{equation}{\|c_1\|} = {\left[ {\int_0^\infty {c_1^2\exp\left( {\frac{{{\zeta^2}}}{4}} \right)\textrm{d}\zeta } } \right]^{1/2}} = \sqrt {{\boldsymbol{a}^H}\boldsymbol{a}} ,\end{equation}

Based on this quantity, one can define the growth rate $\sigma $ as

(3.22)\begin{equation}\sigma = \frac{1}{{{\|c_1\|}}}\frac{{\textrm{d}{\|c_1\|}}}{{\textrm{d}\tau }} = \frac{1}{{2{\boldsymbol{a}^H}\boldsymbol{a}}}\left( {\frac{{\textrm{d}{\boldsymbol{a}^H}}}{{\textrm{d}\tau }}\boldsymbol{a} + {\boldsymbol{a}^H}\frac{{\textrm{d}\boldsymbol{a}}}{{\textrm{d}\tau }}} \right).\end{equation}

From (3.19a) and the basic matrix operation on it, the growth rate defined in (3.22) can be rewritten as

(3.23)\begin{equation}\sigma \tau = \frac{{{\boldsymbol{a}^H}\boldsymbol{Ea}}}{{{\boldsymbol{a}^H}\boldsymbol{a}}},\end{equation}

where $\boldsymbol{E} = ({\boldsymbol{B}^{\boldsymbol{H}}} + \boldsymbol{B})/2$. It should be noted that $\boldsymbol{E} = \boldsymbol{B}$ since the matrix $\boldsymbol{B}$ is normal. For the normal matrix $\boldsymbol{B}$, the following Rayleigh quotient can be defined (Amundson 1966):

(3.24)\begin{equation}R(\boldsymbol{B},\boldsymbol{a}) = \frac{{{\boldsymbol{a}^H}\boldsymbol{Ea}}}{{{\boldsymbol{a}^H}\boldsymbol{a}}}.\end{equation}

The maximum value of $R(\boldsymbol{B},\boldsymbol{a})$ reaches the maximum eigenvalue of $\boldsymbol{B}$ when $\boldsymbol{a}$ is the corresponding eigenvector, i.e.

(3.25)\begin{equation}R(\boldsymbol{B},\boldsymbol{a}) \le \textrm{max\{eig(}\boldsymbol{B}\textrm{)\} }.\end{equation}

By using (3.23) and (3.25), the following relation is obtained:

(3.26)\begin{equation}\sigma \tau \le \textrm{max\{eig(}\boldsymbol{B}\textrm{)\} }.\end{equation}

Therefore, the growth rate of the most dangerous disturbances at a given time $\tau $ can be expressed as

(3.27)\begin{equation}\textrm{max(}\sigma \tau \textrm{)} = \textrm{max\{eig(}\boldsymbol{B}\textrm{)\} }.\end{equation}

For the normal system, where $\boldsymbol{B} = {\boldsymbol{B}^T}$, the detailed procedure to arrive at (3.27) was discussed by Kim & Choi (2015) and Kim (2021). According to the normal mode analysis, for a given $\tau $, the maximum growth rate can be obtained from the maximum eigenvalue of the time-dependent characteristic matrix $\boldsymbol{B}$.

Except for the limiting case of $\chi = 1$, $({C_{sat}}/{\rho _S}) \to 0$ and ${\bar{\mu }_0} = 1$, a fully analytic solution is not possible due to the complex boundary conditions (3.7). However, Kim (2016) proved that the normal mode analysis can be applied to the moving interface problem, where $({C_{sat}}/{\rho _S}) \ne 0$ and the eigenfunctions ${\phi _i}$ and ${\psi _i}$ should be obtained numerically.

3.2.2. For the limiting case of $Da \to 0$

For this reaction-controlled system, the stability equations (3.5)–(3.7) are reduced as

(3.28)\begin{gather}\left( {\frac{{{\partial^2}}}{{\partial {\zeta^2}}} - {k^{{\ast} 2}}} \right)\left\{ {{{\bar{\mu }}_0}\left( {\frac{{{\partial^2}}}{{\partial {\zeta^2}}} - {k^{{\ast} 2}}} \right)w_1^\ast } \right\} = {k^{{\ast} 2}}c_1^\ast ,\end{gather}

(3.29)\begin{gather}\tau \frac{{\partial {c_1}}}{{\partial \tau }} + {w_0}\frac{{\partial {c_1}}}{{\partial \zeta }} + Ra_D^\ast {w_1}\frac{{\textrm{d}({{\tilde{c}}_0}/\sqrt \tau )}}{{\textrm{d}\zeta }} = ({\mathrm{{\mathcal{L}}}_\zeta } - {k^{{\ast} 2}})c_1^\ast ,\end{gather}

under the following boundary conditions:

(3.30a)\begin{gather}w_1^\ast= \frac{{\partial w_1^\ast }}{{\partial \zeta }} = \frac{{\partial c_1^\ast }}{{\partial \zeta }}\;\textrm{at}\;\zeta = 0,\end{gather}

(3.30b)\begin{gather}w_1^\ast \to 0,\quad \frac{{\partial w_1^\ast }}{{\partial \zeta }} \to 0\quad \textrm{and}\quad c_1^\ast \to 0\;\textrm{as}\;\zeta \to \infty ,\end{gather}

where $Ra_D^\ast ( = R{a^\ast }D{a^\ast }) = R{a_D}{\tau ^2}$, $c_1^\ast= {c_1}/Da$ and ${\tilde{c}_0}/\sqrt \tau = 2ierfc(\zeta /2)$.

For the limiting case of ${\bar{\mu }_0} = 1$, the analytic solutions of $c_1^\ast $ and $w_1^\ast $ are given in Appendix B. Even for this limiting case, since we cannot prove the characteristic stability matrix is normal, we should depend on the quasi-steady state approximation (QSSA)

(3.31)\begin{equation}[w_1^\ast ,{c_1}] = [{\bar{w}_1}({k_\ast },\zeta ),{\bar{c}_1}(\zeta )]\textrm{exp(}{\sigma _\ast }{\tau _\ast }\textrm{)}\;\textrm{at}\;\tau = {\tau _\ast },\end{equation}

under the assumption that the temporal growth rate of the base-concentration field is smaller than those of the disturbance quantities. In this equation, ${k_\ast } = k\sqrt {{\tau _\ast }} $. This QSSA implies that all of the time dependencies are frozen at a given time ${\tau _\ast }$, and then $\partial {c_1}/\partial \tau {|_{{\tau _\ast }}} = {\sigma _\ast }{c_1}$ is assumed.

3.2.3. General case

Since the normal mode analysis can be used only for the limiting case of $D{a^\ast } \to \infty $ and $\chi = 1$, the QSSA is applied to solve the stability equations. By applying the QSSA, (3.31), into (3.5)–(3.7), we obtain the following stability equations:

(3.32)\begin{gather}\left\{ {{{\bar{\mu }}_0}\left( {\frac{{{\textrm{d}^2}{{\bar{w}}_1}}}{{\textrm{d}{\zeta^2}}} - k_\ast^2} \right) + 2\frac{{\partial {{\bar{\mu }}_0}}}{{\partial \zeta }}} \right\}\left( {\frac{{{\textrm{d}^2}{w_1}}}{{\textrm{d}{\zeta^2}}} - k_\ast^2} \right){\bar{w}_1} - \frac{{{\partial ^2}{{\bar{\mu }}_0}}}{{\partial {\zeta ^2}}}\left( {\frac{{{\textrm{d}^2}{w_1}}}{{\textrm{d}{\zeta^2}}} + k_\ast^2} \right){\bar{w}_1} = k_\ast ^2{\bar{c}_1},\end{gather}

(3.33)\begin{gather}{\sigma _\ast }\tau {\bar{c}_1} + R{a_\ast }{\bar{w}_1}\frac{{\textrm{d}{c_0}}}{{\textrm{d}\zeta }} = \left\{ {\frac{{{\textrm{d}^2}}}{{\textrm{d}{\zeta^2}}} + \left( {\frac{\zeta }{2} + \frac{{(1 - \chi )\eta }}{2}} \right)\frac{\textrm{d}}{{\textrm{d}\zeta }} - k_\ast^2} \right\}{\bar{c}_1},\end{gather}

under the boundary conditions

(3.34a)\begin{gather}{\bar{w}_1} = \frac{{(1 - \chi )}}{{R{a_\ast }}}D{a_\ast }\left( {\frac{{{C_{sat}}}}{{{\rho_S}}}} \right){\bar{c}_1}\quad \textrm{and}\quad \frac{{\textrm{d}{{\bar{w}}_1}}}{{\textrm{d}\zeta }} = 0\;\textrm{at}\;\zeta =- \eta ,\end{gather}

(3.34b)\begin{gather}\frac{{\textrm{d}{{\bar{c}}_1}}}{{\textrm{d}\zeta }} = D{a_\ast }{\bar{c}_1}\{ 1 - \chi ({C_{sat}}/{\rho _S}){c_0}\} + D{a_\ast }(1 - {c_0})\chi ({C_{sat}}/{\rho _S}){\bar{c}_1}\;\textrm{at}\;\zeta =- \eta ,\end{gather}

(3.34c)\begin{gather}{\bar{w}_1} \to 0,\quad \; \frac{{\textrm{d}{{\bar{w}}_1}}}{{\textrm{d}\zeta }} \to 0\quad \textrm{and}\quad {\bar{c}_1} \to 0\;\textrm{as}\;\zeta \to \infty ,\end{gather}

where $R{a_\ast } = Ra\tau _\ast ^{3/2}$ and $D{a_\ast } = Da\tau _\ast ^{1/2}$. The subscript ‘*’ means the time-dependent quantity which is frozen at $\tau = {\tau _\ast }$, i.e. ${k_\ast }$, $R{a_\ast }$ and $D{a_\ast }$ correspond to ${k^\mathrm{\ast }}$, $R{a^\ast }$ and $D{a^\ast }$ frozen at $\tau = {\tau _\ast }$. The above stability equations (3.32)–(3.34) are to be solved by employing the outward shooting scheme (Kim et al. 2021).

We now will explain how to find the neutral stability curve corresponding to ${\sigma _\ast } = 0$. In order to integrate these stability equations, the proper values of ${\textrm{d}^2}{\bar{w}_1}/\textrm{d}{\zeta ^2}$, ${\textrm{d}^3}{\bar{w}_1}/\textrm{d}{\zeta ^3}$ and ${\bar{c}_1}$ at $\zeta =- \eta $ are assumed for a given $D{a_\ast }$, ${k_\ast }$, $\chi $ and $({C_{sat}}/{\rho _S})$. Since the stability equations and their boundary conditions are all homogeneous, the value of ${\bar{c}_1}$ at $\zeta =- \eta $ can be assigned arbitrarily, and the value of the parameter $R{a_\ast }$ is assumed. This procedure can be understood easily by considering the characteristics of eigenvalue problems. After all of the values at $\zeta =- \eta $ are provided, this eigenvalue problem can be solved numerically. Integration is performed from $\zeta =- \eta $ to an arbitrary lower boundary ${\zeta _f}$ with the fourth-order Runge–Kutta-Gill method. If the guessed value of $R{a_\ast }$, values of ${\textrm{d}^2}{\bar{w}_1}/\textrm{d}{\zeta ^2}$ and ${\textrm{d}^3}{\bar{w}_1}/\textrm{d}{\zeta ^3}$ at $\zeta =- \eta $ are correct, ${\bar{w}_1}\; $, $\textrm{d}{\bar{w}_1}/\textrm{d}\zeta $ and ${\bar{c}_1}$ will vanish at ${\zeta _f}$. To improve the initial guesses, the Newton–Raphson iteration method is used.

For the present system, the position of the lower boundary is infinity. To determine the asymptotic value of $R{a_\ast }$ for ${\zeta _f} \to \infty $, we increase ${\zeta _f}$ by a fixed increment $\Delta \zeta $ systematically as ${\zeta _{f,i + 1}} = {\zeta _{f,i}} + \Delta {\zeta _f}$, and then calculate $R{a_\ast }({\zeta _{f,i + 1}})$. By assuming that $\Delta R{a_{{\ast} ,i + 1}}\{ = R{a_\ast }({\zeta _{f,i + 1}}) - R{a_\ast }({\zeta _{f,i}})\} $ decreases in geometrical sequence, we calculate the common ratio, r, as

(3.35a,b)\begin{equation}r = \frac{{\textrm{d}R}}{{\textrm{d}{\zeta _f}}} = \frac{{{R_n} - {R_{n - 1}}}}{{\Delta {\zeta _f}}}\quad \textrm{and}\quad {R_n} = \frac{{\{ R{a_\ast }({\zeta _{f,n + 1}}) - R{a_\ast }({\zeta _{f,n}})\} }}{{\{ R{a_\ast }({\zeta _{f,n}})\Delta {\zeta _f}\} }}.\end{equation}

Then, $R{a_\ast }(\infty )$, i.e. $R{a_\ast }$ for the infinite boundary, can be obtained as

(3.36)\begin{align} R{a_\ast }(\infty ) & = R{a_\ast }({\zeta _n}) + \sum\limits_{j = 1}^\infty {\{ R{a_\ast }({\zeta _{n + j}}) - R{a_\ast }({\zeta _{n + j - 1}})\} } \notag\\ & \approx R{a_\ast }({\zeta _n}) + \sum\limits_{j = 1}^\infty {\{ R{a_\ast }({\zeta _n}) - R{a_\ast }({\zeta _{n - 1}}){r^j}\} } \notag\\ & = R{a_\ast }({\zeta _n}) + \{ R{a_\ast }({\zeta _n}) - R{a_\ast }({\zeta _{n - 1}})\} \dfrac{r}{{1 - r}}. \end{align}

For the extreme case of $\chi = 1$, $\eta = 0$ and ${\bar{\mu }_0} = 1$, the analytical approximations summarized in Appendices A and B are compared with the numerical solution in figure 3(a). As shown in this figure, the analytic approximate normal mode solutions and the numerical shooting solution based on the QSSA are in good agreement. For the system whose characteristics matrix $\boldsymbol{B}$ is transient but normal, we can freeze the matrix $\boldsymbol{B}$ at a given time (see (3.27)). This means that the present QSSA represents the normal mode solution for this extreme case. This interesting finding cannot be reproduced in the stability equations (3.1)–(3.4) in the $(\tau ,z)$-domain; however, the stability equations (3.5)–(3.7) can be solved in the $(\tau ,\zeta )$-domain. In addition, for $Da \to 0$ ${\bar{\mu }_0} = 1$, the analytically and numerically obtained neutral stability curves based on the QSSA are compared in figure 3(b). As shown in this figure, our numerical solution explains the analytic approximations quite well. Based on this figure, for the general case where the analytical approach is not possible, we obtained the numerical shooting solutions. In addition, the distributions of disturbances are summarized in figures 3(c) and 3(d). The critical conditions of the onset of convection are determined by the minimum point of the neutral stability curves given in figures 3(a) and 3(b). As shown in figures 3(c) and 3(d), the concentration disturbances exist in a narrower region than the velocity disturbances, because for the limiting case of $Sc \gg 1$, the velocity boundary-layer thickness is much larger than the concentration boundary-layer thickness. As discussed previously, because the boundary conditions of base-concentration fields are complex function of $D{a^\ast }$, $({C_{sat}}/{\rho _S})$ and $\chi $, we cannot solve the case field equations analytically. Due to this mathematical difficulty, in the present study, a stability analysis for finite $Da$ case is attempted for the static system of $\chi \to 1$ and ${C_{sat}}/{\rho _s} \to 0$. Instead, the effect of $Da$ on the gravitational instability in a dissolution system will be studied using fully nonlinear numerical simulations.

Figure 3. Neutral stability curves for (a) the transport-controlled system and (b) the reaction-controlled system. The critical conditions given in these figures are determined by the minimum values of the 3-term approximation (a), and 5-term approximation (b). The distributions of disturbances at these critical conditions of (c) the transport-controlled system and (d) the reaction-controlled system are given. Here, an n-term approximation means the analytic solution using n-terms in (3.13), and numerical shooting solutions are obtained using the solution method explained in § 3.2.3. In (c) and (d), disturbance quantities are normalized by the maximum values of ${c_1}$, and red and black colours correspond to the cases in (a) and (b), respectively.

4. Numerical simulations

Since the previous theoretical linear analyses are valid only near the critical time, the convective motion driven by dissolution should be studied by solving (2.12)–(2.17), numerically. In the 2-D domain, we solve (2.12)–(2.17) through an in-house model implementation using COMSOL Multiphysics (2019). The flow field was obtained by using the ‘Laminar Flow’ module of the COMSOL and was used to solve the continuity equation (2.12) and the Navier–Stokes equation (2.13) with the ‘Volume Force’ option for the buoyancy force. The mass transport equation (2.14) was modelled by using the ‘Transport of Dilute Species’ module, where the flow fields were adopted from the Navier–Stokes equation. Using the initial condition (2.15a), numerical integrations proceeded. The boundary conditions of (2.15b) are implemented by a leaking wall condition, and (2.15c) is the fixed concentration condition. To implement the salt–liquid interface movement, the ‘Deformed Geometry’ module was used with the prescribed normal velocity given in (2.16). By introducing the normal velocity, (2.16), at the upper boundary through translational wall movement the interface movement was solved simultaneously with the velocity and concentration field. Besides the boundary conditions (2.15), the following periodic conditions were imposed:

(4.1a,b)\begin{equation}p{|_{x =- L/2}} = p{|_{x = L/2}}\quad \textrm{and}\quad c{|_{x =- L/2}} = c{|_{x = L/2}},\end{equation}

where L corresponds to the aspect ratio. Near the dissolving surface, where the gradients are high, we used the ‘boundary layer’ refinement option in constructing a mesh structure. Geometry and boundary conditions are schematized in figure 4.

Figure 4. Schematic diagram of the system geometry and boundary conditions used in this numerical simulation study.

The COMSOL Multiphysics software has a wide variety of options in choosing the time-dependent solver. In the present study, a first- or second-order variable step size backward differentiation formula (BDF) was used. The order of the BDF solver is determined by the degree of the interpolating polynomial. At each time step, the system of nonlinear algebraic equations is linearized employing the Newton method, and the resulting linearized system is solved by the PARDISO direct solver which is fast, robust and multi-core capable. The scaled absolute tolerance factors of 0.05 and 1 and were set for the concentration and the pressure, respectively. By changing the integration time step, the space discretization and the aspect ratio, $AR( = L)$, we tested the convergence of solutions. Further details on the calculations were discussed in Appendix C.

In the present study, we set $L = 1$, i.e. the aspect ratio is 1. We traced the temporal evolution of the concentration field by using the following dissolved quantity and Sherwood numbers:

(4.2a)\begin{gather}D = \int_\varOmega {c\,\textrm{d}\varOmega } = \int_{ - L/2}^{L/2} {\int_h^1 {c\,\textrm{d}z\,\textrm{d}\kern0.06em x} } ,\end{gather}

(4.2b,c)\begin{gather}S{h_T} = \frac{1}{L}\int_{ - L/2}^{L/2} {{{\left. {\frac{{\partial c}}{{\partial z}}} \right|}_h}\textrm{d}\kern0.06em x} \quad \textrm{and}\quad S{h_R} = \frac{1}{L}\int_{ - L/2}^{L/2} {{{\left. {\frac{1}{{\tilde{c}}}} \right|}_h}\textrm{d}\kern0.06em x} ,\end{gather}

where the subscripts ‘T’ and ‘R’ mean the transport-controlled quantity and the reaction-controlled one, respectively. Generally, the Sherwood number, which has the same meaning as the Nusselt number of the heat transfer system, is defined as

(4.3)\begin{equation}S{h_T} = \frac{{\textrm{total}\;\textrm{mass}\;\textrm{transfer}\;\textrm{rate}}}{{\textrm{diffusional}\;\textrm{mass}\;\textrm{transfer}\;\textrm{rate}}} = \frac{{\langle - \mathrm{{\mathcal{D}}}\partial C/\partial Z{|_H}\rangle }}{{\mathrm{{\mathcal{D}}}\Delta C/d}} = \left\langle {{{\left. { - \frac{{\partial c}}{{\partial z}}} \right|}_h}} \right\rangle ,\end{equation}

for the transport-controlled system where the interface concentration is constant. Here, $\langle \cdot \rangle $ means the surface average. However, for the reaction-controlled system, the above Sherwood number should be modified as

(4.4)\begin{equation}S{h_R} = \frac{{ - \mathrm{{\mathcal{D}}}\partial C/\partial Z{|_H}}}{{\langle \mathrm{{\mathcal{D}}}\Delta C/d\rangle }} = \frac{{ - \partial c/\partial z{|_h}}}{{\langle c{|_h} - c{|_{z = 1}}\rangle }} = \frac{{ - \partial c/\partial z{|_h}}}{{\langle c{|_h}\rangle }} = \frac{1}{{\langle \tilde{c}{|_0}\rangle }},\end{equation}

where $\tilde{c} = c/Da$ because $h \to 0$ and $\partial c/\partial z{|_0} =- Da$ are given by the boundary condition (2.31).

For the example case of $Sc = {10^3}$ and fixed boundary system, i.e. ${C_{sat}} \to 0$ and $\chi = 1$, the above quantities are summarized in figures 5 and 6. For the limiting cases of $Da \to \infty $ and $Da \to 0$, by combining equations (2.22) (2.35) and (4.2) we can get the following results:

(4.5a)\begin{gather}{D_{0,T}} = \int_\varOmega {{c_0}\,\textrm{d}\varOmega } = \sqrt \tau \int_{ - L/2}^{L/2} {\int_h^1 {{c_0}\,\textrm{d}\zeta \,\textrm{d}\kern0.06em x} } = \frac{2}{{erfc( - \chi \eta /2)}}ierfc( - \chi \eta /2)L\sqrt \tau ,\end{gather}

(4.5b)\begin{gather}S{h_{T,0}} = \frac{1}{L}\int_{ - L/2}^{L/2} {{{\left. {\frac{{\partial {c_0}}}{{\partial z}}} \right|}_h}\textrm{d}\kern0.06em x} = \frac{1}{{\sqrt {{\rm \pi} \tau } }}\frac{{\exp ( - {\eta ^2}{\chi ^2}/4)}}{{erfc( - \eta \chi /2)}},\end{gather}

(4.5c)\begin{gather}S{h_{R,0}} = \frac{1}{L}\int_{ - L/2}^{L/2} {{{\left. {\frac{1}{{\tilde{c}}}} \right|}_h}\textrm{d}\kern0.06em x} = \frac{1}{{ierfc(0)}}\frac{1}{{\sqrt {4\tau } }}.\end{gather}

As shown in figures 5 and 6, the present analytic and numerical solutions are in good agreement until the Sherwood numbers defined in (4.2b) and (4.2c) reach their minimum values. Moreover, we identified the onset time of instability as the time ${\tau _d}$ at which the Sherwood numbers show their minimum. Figures 5 and 6 show that nonlinear dynamics, such as the mass transfer enhancement and the pattern formation on the dissolving surface, is not observed during not only the pure diffusion region, $0 < \tau < {\tau _c}$, but also the diffusion-dominant region, ${\tau _c} < \tau < {\tau _d}$. In the present numerical simulations, we will focus on the effect of the various physicochemical parameters on the nonlinear dynamics for the convection-dominant region, $\tau > {\tau _d}$.

Figure 5. Numerical simulation results for the transport-controlled fixed boundary system: (a) temporal evolution of the dissolved quantity $(Ra^{1/3}D)$ and (b) temporal evolution of the mass transfer rate $ (Ra^{-1/3}Sh_T)$. For the case of $Ra \ge {10^5}$, the deviation time ${\tau _d}$, from which the mass transfer flux starts to deviate from the diffusional flux, can be scaled as $R{a^{2/3}}{\tau _d}\sim \textrm{const}$.

Figure 6. Numerical simulation results for the reaction-controlled system. For the case of $R{a_D} \ge {10^5}$, the deviation time ${\tau _d}$, from which the mass transfer flux $( Ra_{D}^{-1/4} Sh_{R})$ starts to deviate from the diffusional flux, can be scaled as $Ra_D^{1/2}{\tau _d}\sim \textrm{const}$.

5. Results and discussion

5.1. Asymptotic scaling relationships

At this point, we study the effect of the Rayleigh number on the onset and growth of dissolution-driven flow instability. For a fixed interface system, from the critical conditions $Ra_c^\ast ( = Ra\tau _c^{3/2}) = 20.73$ and $k_c^\ast ( = k\tau _c^{1/2}) = 0.53$ (see figure 3a), and $Ra_{D,c}^\ast ( = R{a_D}\tau _c^2) = 20.07$ and $k_c^\ast ( = k\tau _c^{1/2}) = 0.53$ (see figure 3b), the following critical conditions can be obtained from a linear stability analysis:

(5.1a)\begin{gather}{\tau _c} = 7.54R{a^{ - 2/3}}\quad \textrm{and}\quad {k_c}( = 2{\rm \pi} d/{\lambda _c}) = 0.197R{a^{1/3}}\notag\\ \textrm{for}\;Da \to \infty ,\;Ra \to \infty \;\textrm{and}\;\varGamma = 0,\end{gather}

(5.1b)\begin{gather}{\tau _c} = 4.48Ra_D^{ - 1/2}\quad \textrm{and}\quad {k_c}( = 2{\rm \pi} d/{\lambda _c}) = 0.250Ra_D^{1/4}\; \notag\\ \textrm{for}\;Da \to 0,\;DaRa \to \infty \;\textrm{and}\;\varGamma = 0,\end{gather}

where ${\lambda _c}$ is the critical wavelength. In addition, from the numerical simulations in figures 5 and 6, there exists a deviation time, ${\tau _d}$, from which mass transfer characteristics start to deviate from those of pure diffusion states. In addition, figures 5 and 6 show that the mass transfer characteristics deviate from those of pure diffusion state at a certain time which is indicated as a deviation time,$\; {\tau _d}$. In the present study, the deviation time, ${\tau _d}$ is defined as the point at which the Sherwood number shows its first local minimum. Even though there exists a certain factor between ${\tau _c}$ from the linear theory and ${\tau _d}$ form the numerical simulations, the deviation time, ${\tau _d}$, follows scaling relationships identical to those from the linear stability analysis, i.e. ${\tau _d}\sim R{a^{ - 2/3}}$ when $Da \to \infty $ and $Ra > {10^5}$ (see figure 5) and ${\tau _d}\sim Ra_D^{1/4}$ when $Da \to 0$ and $R{a_D} > {10^5}$ (see figure 6). The difference between ${\tau _c}$ and ${\tau _d}$ will be discussed further in § 5.4. For the sake of consistency, we used the initial depth of the solution, d, as the length scale in the present linear stability analyses and numerical simulations.

For a salt-dissolving system $(Da \gg 1)$, Philippi et al. (2019) and Berhanu et al. (2021) proposed scaling relationships identical to those here, i.e. (5.1). They suggested the following relationships:

(5.2a,b)\begin{equation}{\tau _c}\sim R{a^{ - 2/3}}\quad \textrm{and}\quad {\lambda _c}\sim R{a^{ - 1/3}}.\end{equation}

Philippi et al. (2019) also conducted experiments and numerical simulations. Based on their salt-dissolution experiments, the emission time of the first plume $({\tau _{em}})$ is correlated as

(5.3)\begin{equation}{\tau _{em}} = 12.8R{a^{ - 2/3}}.\end{equation}

These scaling relationships are identical to the linear and nonlinear analyses in this study. For the other limiting case of $Da \ll 1$, Philippi et al. (2019) also showed the following scaling relationships:

(5.4a,b)\begin{equation}{\tau _c}\sim Ra_D^{ - 1/2}\quad \textrm{and}\quad {\lambda _c}\sim Ra_D^{ - 1/4},\end{equation}

which are identical to (5.1b).

5.2. Effect of Damköhler number

The effect of the Damköhler number Da on the dissolution-driven instability is studied. Depending on Da, the dissolution behavior shows a transition between the transport-controlled regime $(Da \gg 1)$ and the reaction-controlled regime $(Da \ll 1)$.

For the example case of $Sc = {10^3}$, $\chi = 0.75$, $\mathrm{\Delta }C/{\rho _S} = 0.1$, $\varGamma = 0$ and $Ra = {10^6}$, the effects of $Da$ on the mass transfer characteristics and the concentration field are summarized in figures 7 and 8, respectively. As shown in figure 7, in the transport-controlled regime $(Da \gg 1)$, the surface average concentration on the dissolving surface, $\langle {c_i}\rangle $, is strongly dependent on $Da$ and $\tau $ (see figure 7a). However, because for the cases of $Da \ge {10^3}$, ${c_s}$ deviates from 1 only for the small $\tau $-region, the dissolved quantity of salt, D, increases with an increase in $Da$ and reaches saturation, and the mass transfer rate, $\textrm{S}{\textrm{h}_T}$, is nearly equal to that for the case of $Da \to \infty $. Based on these findings, for the system of $Da \ge {10^3}$, the nonlinear dynamics of dissolution can be approximated by that of the transport-controlled system. For the limiting case of $Da \to \infty $, before the onset of instability, the present numerical simulation (black solid lines in figure 7b,c) is in good agreement with its analytic value (black dashed lines in figure 7b,c). From the snapshots of the concentration fields given in figure 7(d), it seems that, with increasing $Da$, the wavelength of the instability decreases i.e. more downward plums are expected. In the reaction-controlled regime $(Da \ll 1)$, from figure 8(a), the proper governing parameter is $R{a_D}( = RaDa)$ rather than $Ra$ or $Da$, and from figure 8(b), any critical difference in the mass transfer rate, $S{h_R}$, and the concentration distribution are not observed for the case of $Da < 0.1$.

Figure 7. Effects of Da on (a) the temporal evolution of surface concentration $(c_s)$ for the high-Da system, (b) the temporal evolution outcomes of the dissolved quantity $(Ra^{1/3}D)$ for the high-Da system, (c) the temporal evolutions of the mass transfer rate for the high-Da system and (d) the snapshots of the concentration field at $R{a^{2/3}}\tau = 100$. For the present system, the value of $Da$ plays little role in the nonlinear dynamics of dissolution and, therefore, the salt concentration on the dissolving surface, ${C_i}$, can be approximated by ${C_{sat}}$ if $Da \ge {10^3}$.

Figure 8. Effects of Da on (a) the temporal evolution of the mass transfer rate $(Ra_{D}^{-1/4}Sh)$ for the low-Da system and (b) the snapshots of the concentration field at $Ra_D^{1/2}\tau = 100$. For the slow reaction system of $Da \le 0.1$, the nonlinear dynamics of the dissolution is governed by $R{a_D}( = RaDa)$ rather than $Ra$ and $Da$ and no critical difference can be observed among these systems.

Even though Berhanu et al. (2021) conducted linear stability analysis for the large $Da$ system, a systematic stability analysis to reveal the effect of $Da$ on the onset of instability has not been tried. For Philippi et al.'s (2019) and Berhanu et al.'s (2021) static system, the effect of $Da$ on the onset of instability is summarized in figure 9. For the high $Da$ system, the deviation time, ${\tau _d}$, as the time at which $S{h_T}$ (see (3.36b)) its minimum is compared with ${\tau _c}$ of the present linear stability analysis. In addition, for the low $Da$ system, ${\tau _d}$ corresponding to the minimum point of $S{h_R}\; $ is compared with ${\tau _c}$ in terms of $R{a_D}( = RaDa)$. Figure 9 clearly shows that the present theoretical predictions and numerical simulations follow scaling relationships based on the linear stability analysis, i.e. (5.1). As summarized in figure 9, the present theoretical and numerical studies show that, for the intermediate Da system, $1 \le Da \le 100$, the Damköhler number plays an important role in the dissolution-driven convective instability, and therefore, both the diffusional mass transfer and the surface dissolution reaction should be considered. This finding is quite similar to Colombani & Bert's (2007) proposal that, for the mixed kinetics regime where $0.1 \le Da( = {\alpha ^{ - 1}}) \le 100$, the dissolution of crystal is controlled by a mixed surface reaction and diffusion kinetics. It has been known that there exists a certain growth period during which the most unstable instability onset at ${\tau _c}$ grows enough to be detected experimentally (Kim & Choi 2006). As shown in figure 9, there is the consistent relationship between the numerically determined ${\tau _d}$ and the theoretically obtained ${\tau _c}$, i.e. ${\tau _d} - {\tau _c} = 5{\tau _c}$.

Figure 9. Effect of $Da$ on the deviation time ${\tau _d}$ for the case of $Sc = {10^3}$, $\chi \to 1$, $\Delta C/{\rho _S} \to 0$, $\varGamma = 0$ and $Ra = {10^6}$. Theoretical and numerical results are plotted in the (a) high and (b) low ranges of $Da$. The parameter ${\tau _c}$ is the onset time determined from the linear stability analysis and ${\tau _d}$ is the deviation time when the mass transfer rate deviates from its diffusional limit. Both the linear stability analysis and the numerical simulation outcomes are qualitatively in good agreement.

5.3. Effect of $\chi $ and ${C_{sat}}/{\rho _S}$

As discussed in § 2.1, during the diffusion-only period, $\chi $ accelerates the dissolution of the salt for the transport-controlled system. The effect of $\chi $ on the dissolution rate is important, especially for a high ${C_{sat}}/{\rho _S}$-system. However, for the reaction-controlled system, i.e. $Da \to 0$, both $\chi $ and ${C_{sat}}/{\rho _S}\; $ barely affect the onset of instability and the dissolution rate because, in this limiting case, the most dominant parameter is $RaDa$. Herein, the effects of $\chi $ and ${C_{sat}}/{\rho _S}$ on the onset of instability and the dissolution rate are investigated while focusing on the transport-controlled system. For the case of $Da \to \infty $ and ${C_{sat}}/{\rho _S} = 0.3$, the effect of $\chi $ on the enhancement of the dissolution rate is summarized in figure 10. As shown in figure 10(a), $\chi $ slightly retards the mass transfer of solute at the dissolving surface. However, the effects of $\chi $ on the dissolution quantity D, the change of the heigh $\Delta {h_{avg}}$ and the concentration fields are not critical as time goes on, as shown in figure 10(b–d). The effects of $\chi $ on the characteristic times are compared in figure 11. The critical times predicted by the linear and nonlinear analyses show a similar tendency on $\chi $; i.e. $\chi $ makes the system unstable slightly. In addition, for the case of $Da \to \infty $ and $\chi = 1$, the effects of ${C_{sat}}/{\rho _S}$ on the dissolution rate, the enhanced mass transfer by the instability and the concentration fields are summarized in figure 12. It is clear that ${C_{sat}}/{\rho _S}$ makes the system unstable and accelerates the onset of instability motion. The effect of ${C_{sat}}/{\rho _S}$ on the critical time predicted by the linear and nonlinear analyses is summarized in figure 13. Both the linear and the nonlinear analyses show the same tendency, i.e. ${C_{sat}}/{\rho _S}$ makes the system slightly unstable. As discussed in figure 9, it seems that there exists a certain growth period, $5{\tau _c}$, and during this certain time period, the most unstable instability onset at ${\tau _c}$ grows enough to enhance the mass transfer rate (see figures 11 and 13).

Figure 10. Effect of $\chi $ on (a) temporal evolution of the mass transfer rate $(Ra^{-1/3}Sh_T)$, (b) the dissolved quantity $({Ra^{1/3}}D)$, (c) the change of height $({Ra^{1/3}}\Delta h_{avg})$ and (d) concentration field at $R{a^{2/3}}\tau = 80$ for the transport-controlled system. Dashed lines correspond to the diffusional solutions summarized in (3.32).

Figure 11. Comparison of (a) the critical times $({Ra^{2/3}}\tau)$ obtained from the linear stability analysis and the numerical simulation at various $\chi $-values, and (b) the critical wavenumber $(k/Ra^{1/3})$ with Sullivan et al.'s (1996) experiments. Present theoretical predictions explain the numerical simulations and the previous experiments qualitatively.

Figure 12. Effect of ${C_{sat}}/{\rho _S}$ on (a) the dissolved quantity $(Ra^{1/3}D)$, (b) the mass transfer rate $(Ra^{-1/3}Sh_T)$ and (c) concentration profile at $R{a^{2/3}}\tau = 50$ for the transport-controlled system, i.e. $Da \to \infty $. Dashed lines correspond to the diffusional solutions summarized in (3.32).

Figure 13. Comparison of (a) the critical times $({Ra^{2/3}}\tau)$ obtained from the linear stability analysis and the numerical simulation at various ${C_{sat}}/{\rho _S}$-values, and (b) the critical wavenumber $(k/Ra^{2/3})$ with Sullivan et al.'s (1996) experimental results for NaCl (, ${C_{sat}}/{\rho _S} = 0.135$), KBr ($\Delta$, ${C_{sat}}/{\rho _S} = 0.247$) and KCl (▾, ${C_{sat}}/{\rho _S} = 0.171$). Present theoretical predictions explain the numerical simulations and the previous experiments qualitatively.

Sullivan et al. (1996) experimentally studied pattern formation on the dissolving surface of salts, in their case NaCl(${C_{sat}}/{\rho _S} = 0.135$), KBr(${C_{sat}}/{\rho _S} = 0.247$) and KCl(${C_{sat}}/{\rho _S} = 0.171$). According to their experimental visualization results, the patterns form on the dissolving interface differently depending on the dissolving species. In Sullivan et al.'s (1996) dissolution experiments, they suggested that the dimensionless wavenumber of the dissolution pattern has the following relation:

(5.5)\begin{equation}{k_c} \approx 0.065R{a^{1/3}},\end{equation}

regardless of the dissolving salt. This means that, in real experimental conditions, the Rayleigh number is the most important parameter to control the wavelength of the dissolving pattern.

In figures 11 and 13, the effects of $\chi $ and ${C_{sat}}/{\rho _S}$ on the critical wavenumber are summarized and compared with Sullivan et al.'s (1996) experimental study. Even though the present linear analysis reproduced Sullivan et al.'s (1996) experiments qualitatively, the predicted wavenumber is nearly three times larger than the experimental results. In other words, the experimentally observed wavelengths are three times smaller than the present prediction. Generally, the linear theory predicted a critical wavenumber in good agreement with experimental outcome if it was measured near the onset condition (Choi, Park & Kim 2004; Kim, Park & Choi 2005; Kim, Choi & Yeo 2007). However, in Sullivan et al.'s study (1996), they measured the wavenumber at the quasi-stationary regime far from the onset, explaining that the discrepancy between their experiments and linear stability analysis is attributed to the difference between experimental measuring time and theoretical onset condition. Experimentally determined wavenumbers are insensitive to ${C_{sat}}/{\rho _S}$. The same trend can be expected from the present linear analysis (see figure 13b).

Unlike the present study, Philippi et al. (2019) conducted a numerical simulation under the following boundary conditions:

(5.6)\begin{equation}\boldsymbol{u} = 0\quad \textrm{and}\quad \boldsymbol{\nabla }c\boldsymbol{\cdot }\boldsymbol{n} =- Da(1 - c)\;\textrm{at}\;z = h = 0,\end{equation}

which correspond to the limiting case of $\chi = 1$ and ${C_{sat}}/{\rho _S} = 0$. Through this static assumption, the dissolving surface is fixed and pattern formation on the dissolving surface is therefore not allowed. Therefore, by comparing this study that introduces moving interface conditions, (2.15)–(2.17), with a previous study that uses a fixed interface under the static assumption, (5.6), the influence of the interface condition on the pattern formation outcome can be distinguished.

5.4. Effect of viscosity variation

For the thermal system, Davaille & Jaupart (1993, 1994) experimentally studied the effect of thermally induced viscosity decrease on the onset time of instability, which was theoretically explained by Korenaga & Jordan (2003, 2004) and Kim & Choi (2006) by considering temperature-dependent viscosity variations. Through these previous experimental and theoretical studies, it is known that the thermally induced viscosity decrease accelerates the onset of instability motion. Here, we study the effect of the concentration-dependent viscosity thickening on the onset and the growth of instability motions.

For various $\varGamma $-values, numerical predictions of the dissolved quantity and mass transfer rate are summarized in figure 14. As discussed in previous studies (Davaille & Jaupart 1993; Korenaga & Jordan 2003; Kim & Choi 2006; Sabet et al. 2021; Kim 2022), figure 14 shows that the thickening of the viscosity impedes the onset and growth of instability motions. The critical time obtained from the linear stability analysis is also compared with the present numerical simulation results in figure 15. As discussed in figures 9, 11 and 13, there exists a certain difference between ${\tau _c}$, which is predicted by the linear stability theory, and ${\tau _d}$, which is the deviation time determined through the nonlinear numerical simulations. A similar trend also exists between ${\tau _c}$ and the experimentally determined characteristic time, such as Philippi et al.'s (2019) ${\tau _{em}}$, at which the first plume is emitted. These differences between the linear theory and numerical simulation or experiment have been understood as the growth period for the infinitesimal disturbance growing exponentially. During this growth period, the diffusive forces of mass diffusion and viscosity suppress the break of the boundary layer. Then, beyond this growth period, the nonlinear dynamics of instabilities can be observed experimentally or via numerical simulations.

Figure 14. Effect of $\varGamma $ on (a) the dissolved quantity $(DRa^{1/3})$ and (b) the mass transfer rate $({Sh_T} /Ra^{1/3})$.

Figure 15. Comparison of the critical times $(\tau Ra^{2/3})$ obtained from the linear stability analysis and the numerical simulation at various $\varGamma $-values.

Unlike the constant viscosity system, where the growth period is expected to be $5{\tau _c}$, the difference between the linear stability theory and the numerical simulation increases with an increase in the $\varGamma $-value; i.e. the viscosity-thickening effect increases. As shown in figure 5(b), for the case of the low $Ra$, the deviation time, ${\tau _d}$, does not follow the scaling relation for the high $Ra$ system, because the penetration depth at ${t_d}$, $\Delta ({t_d})$, is comparable to the depth of the liquid layer, d, and then, the deep-pool approximation no longer holds. Thus, for cases with a high $\varGamma $-value, the onset of instability is retarded and the bottom boundary effect can therefore further impede the onset of motion in the numerical simulation, whereas this bottom boundary effect is not considered in the linear stability analysis. If the $\varGamma $-value itself determines the growth period, the relation ${\tau _d} - {\tau _c} = 5{\tau _c}$ should not hold for the case of $\varGamma \ne 0$. However, the relation ${\tau _d} - {\tau _c} = 5{\tau _c}$ does hold even for the viscosity-thinning system, i.e. $\varGamma =- 1$. This implies that the important factor in determining the growth period is the relative thickness of the penetration depth rather than the viscosity-thickening effect, and the relation ${\tau _d} - {\tau _c} = 5{\tau _c}$ holds for the deep-pool system, where the dimensionless penetration depth at ${\tau _d}$, $\delta ({\tau _d}) \ll 1$.

5.5. Pattern formation on the dissolving surface

As discussed in § 5.3, the static assumption of Philippi et al. (2019), (5.1), cannot allow the movement of a salt-solution interface. Hence, we cannot computationally visualize the pattern formation at the solid salt–liquid interface as observed experimentally by Sullivan et al. (1996). To investigate the formation of the dissolution pattern, the dissolution of the NaCl horizontal salt into water is taken as a model system in this study, with the corresponding physical properties employed as shown in table 1 for the 3-D numerical simulations. Because the dissolution rate constant of a NaCl–water system is large, we assume that $Da$ of the present system is large, i.e. the present system is a transport-controlled one. For $Ra = {10^7}$, the temporal evolution of the height of the liquid layer is summarized in figure 16, where the temporal evolution of the dissolving pattern is clearly shown after a purely diffusion-dominant period (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.85). This pattern formation cannot be observed under the static assumption where $\partial h/\partial \tau = 0$ is assumed. In addition, the effects of the Rayleigh number on the concentration field and the pattern on the dissolving surface are summarized in figure 17. As predicted by the linear stability analysis, the wavelength of the instability motion decreases with an increase in the Rayleigh number, a result identical to the experimental finding of Sullivan et al. (1996). Because the numerical simulations in the present work was conducted with enough patterns (see figures 16 and 17), we can observe that the influence of the lateral boundary condition on the development of the instability pattern becomes less important with increasing $Ra$. Furthermore, as summarized in figure 16 and the supplementary materials, more instability patterns are developed, i.e. the wavenumber of the instability pattern increases as time goes on. Even though we cannot extract the dominant wavenumber from the numerical simulations, the present numerical simulations may relax the discrepancy in wavenumber shown in figures 11 and 13.

Table 1. List of symbols used in this study and their values for the model system of NaCl–salt and the solution. All data are from Philippi et al. (2019) and references therein.

Figure 16. Temporal evolution of the pattern on the dissolving interface of NaCl–water system at $Ra = {10^7}$, $Da \to \infty $, ${c_{sat}}/{\rho _S} = 0.135$ and $\varGamma = 0.688$. The movie for this case available as supplementary movie 1.

Figure 17. Effect of the Rayleigh number on the concentration field (upper panels) and the pattern formed on the dissolving surface (lower panels) at $\tau R{a^{2/3}} = 300$, for the case of $Da \to \infty $, ${c_{sat}}/{\rho _S} = 0.135$ and $\varGamma = 0.688$.

From the interfacial conditions (2.16a) and (2.16b), the dissolution rate can be determined by the concentration distribution at the interface if the parameters, $Da$, $\chi $ and ${C_{sat}}/{\rho _S}$ are given. Furthermore, the gravitational instability induces the non-uniformity of the concentration field, resulting in an inhomogeneous dissolution rate. The inhomogeneity of the dissolution rate carves an irregular pattern on the dissolving surface. We can deduce that the Rayleigh number is one of the key factors which control the pattern formation on the dissolving surface. To determine the effect of the Rayleigh number on the pattern formation process, the temporal evolution outcomes of the dissolution rate; D, $\Delta {h_{avg}}$ and ${\sigma _h}$, are correspondingly summarized in figures 18 and 19. Here, $\Delta {h_{avg}}$ and ${\sigma _h}$ are the average value of the change in the height of the dissolving interface and its standard deviation, respectively. As expected from (4.5a), figure 18 predicts the relationship of $D\sim \sqrt \tau $ for the diffusion-dominant regime. As shown in figure 19, pattern formation driven by salt dissolution evolves depending on the Rayleigh number as follows:

(5.7a)\begin{gather}\Delta {h_{avg}}\sim \sqrt \tau \;\textrm{for}\;\textrm{the}\;\textrm{diffusion-dominant}\;\textrm{regime}\;\textrm{of}\;R{a^{2/3}}\tau \le 40,\end{gather}

(5.7b)\begin{gather}\Delta {h_{avg}}\sim R{a^{0.79/3}}{\tau ^{0.86}}\quad \textrm{and}\quad {\sigma _h}\sim R{a^{0.58/3}}{\tau ^{0.79}}\;\textrm{for}\;\textrm{the}\;\textrm{convection}\;\textrm{regime}\;\textrm{of}\;\notag\\ 40 \le R{a^{2/3}}\tau .\end{gather}

It should be noted that, because we are focused on the initial stage of pattern formation, we did not try to find the upper limit of the convection regime. From figures 18(b) and 19(a), the relations $R{a^{1/3}}D\sim {(R{a^{1/3}}{\tau ^{1/2}})^{1.69}}$ and $R{a^{1/3}}\Delta {h_{avg}}\sim {(R{a^{2/3}}\tau )^{0.86}}$ support each other. If the formation of the pattern on the dissolving surface can be described by the variation of the height of the dissolving interface, the relationship between them will be ${\sigma _h}\sim R{a^{0.58/3}}{\tau ^{0.79}}$ for the initial stage of the pattern formation. For the fully developed turbulent state or the quasi-steady state, the mass transfer flux is constant, suggested by Philippi et al. (2019) as follows:

(5.8)\begin{equation}{V_d}\left( { \approx \frac{{\textrm{d}(\Delta {H_{avg}})}}{{\textrm{d}t}}} \right)\sim {V_{Pattern}}\left( { \approx \frac{{\textrm{d}{\sigma_H}}}{{\textrm{d}t}}} \right)\sim \textrm{const}\textrm{.}\quad \textrm{for}\;Da \gg 1.\end{equation}

This corresponds to $\Delta {h_{avg}}\sim R{a^{1/3}}\tau $ and ${\sigma _h}\sim R{a^{1/3}}\tau $. Furthermore, for the quasi-steady state regime, the relation $\Delta {h_{avg}}\sim R{a^{1/3}}\tau $ has been observed experimentally (Sullivan et al. 1996; Davis Wykes et al. 2018; Cohen et al. 2020). From (5.7) and (5.8) and figure 19, this study proposes the following relations:

(5.9a,b)\begin{equation}R{a^{1/3}}\Delta {h_{avg}}\sim {(R{a^{2/3}}\tau )^m}\quad \textrm{and}\quad R{a^{1/3}}{\sigma _h}\sim {(R{a^{2/3}}\tau )^n}.\end{equation}

The exponents m and n are given in table 2 for each regime. It is interesting that the present numerical simulation results for the convection regime, $m \approx 0.89$ and $n \approx 0.79$, roughly follow Philippi et al.'s (2019) suggestion for the quasi-steady state regime, i.e. $m = n = 1$. From this finding, the present study for the diffusion-dominant and convection regimes complements Philippi et al.'s (2019) work for the quasi-steady state regime. Furthermore, from (5.9), it is found that $d/R{a^{1/3}}$ and ${d^2}/(\mathrm{{\mathcal{D}}}R{a^{2/3}})$ are the appropriate time and length scales in the early stage of pattern formation. Furthermore, we can expect the transition from the diffusion-dominant region to the convection-dominant regime around $R{a^{2/3}}{\tau _t} \approx 40$. As expected, this transition time is closely related to the present theoretical onset time and numerical deviation time, i.e. ${\tau _t} \approx {\tau _d} \approx 6{\tau _c}$.

Figure 18. Effect of the Rayleigh number on the dissolved quantity of NaCl in water: (a) $DR{a^{1/3}}$ vs. $\tau R{a^{1/3}}$ and (b) $DR{a^{1/3}}$ vs. ${\tau ^{1/2}}R{a^{1/3}}$ for the case of $Da \to \infty $, ${c_{sat}}/{\rho _S} = 0.135$ and $\varGamma = 0.688$.

Figure 19. Temporal evolution of (a) the average height of salt $(Ra^{1/3}\Delta h_{avg})$ and (b) the standard deviation of height of salt $(Ra^{1/3}\sigma_h)$ for the case of $Da \to \infty $, ${c_{sat}}/{\rho _S} = 0.135$ and $\varGamma = 0.688$.

Table 2. Exponents m and n for each regime.

^aPredictions from Philippi et al. (2019).

As discussed in § 5.3, the static assumption of Philippi et al. (2019) corresponds to the extreme case of the present system. In the transport-controlled limit, by forcing $\chi = 1$ and ${C_{sat}}/{\rho _S} = 0$, we repeated the 3-D simulations for a NaCl–water system. For the case of $Da \to \infty$ and $Ra = {10^7}$, the temporal evolution of the 3-D concentration fields without the static assumption (the present study) and with the static assumption ($\chi = 1$ and ${C_{sat}}/{\rho _S} = 0$) are given in the supplementary movies 2 and 3, respectively. As expected, any inhomogeneity of the surface concentration cannot be observed in both cases under the condition that the dimensionless concentration of salt on the dissolving surface is $c \to 1$ as $Da \to \infty $ (see (5.6)). This means that the pattern formed on the dissolving surface is due to the inhomogeneity of the mass transfer rate rather than that of the surface concentration. The inhomogeneity of the mass transfer rate is governed by the solute concentration gradient at the dissolving surface or the thickness of concentration boundary layer, because $\boldsymbol{\nabla }C{|_H}\sim \Delta C/{\delta _C}$ where $\boldsymbol{\nabla }C{|_H}$ is the solute concentration gradient at the dissolving surface, $\Delta C$ is the concentration difference between the dissolving surface and the bulk solution and ${\delta _C}$ is the concentration boundary-layer thickness. According to this relation, the boundary-layer thickness is inversely proportional to the concentration gradient at the dissolving surface, i.e. $\boldsymbol{\nabla }C{|_H}\sim 1/{\delta _C}$, because the concentration difference $\Delta C$ is constant in the transport-controlled limit. Because Philippi et al. (2019) explained the pattern formation on the dissolving surface by the inhomogeneity of the surface concentration under the static assumption, their analysis cannot be applicable to the transport-controlled limit. However, the application of the moving interface conditions coupled with dissolution reaction, (2.15b), (2.16a) and (2.16b), clearly shows the development of pattern on the dissolving surface (see supplementary movies 1 and 2). In the present study, the pattern formation on the dissolving salt system is described by using the dimensionless time, the Damköhler number and the Rayleigh number. To compare the present numerical simulation with the real system (NaCl–water system), the Damköhler number, the Rayleigh number and the conversion factor for the dimensionless time to the real time are summarized in table 3 based on the physical properties of the NaCl–water system. Because the Damköhler number and Rayleigh number corresponding to $d = 1\;\textrm{cm}$ are $Da = 3.125 \times {10^3}$ and $Ra = 1.217 \times {10^9}$, respectively, we can use the critical conditions for the case of $\varGamma = 0.688$, $Da \to \infty $ and $d \ge 1\;\textrm{cm}$ (see figure 15). For this practical case, we can predict the following critical and deviation conditions:

(5.10a–c)\begin{equation}{t_c} = 0.50\;\textrm{s,}\quad {t_d} = 3.00\;\textrm{s}\quad \textrm{and}\quad {\lambda _c} = 0.32\;\textrm{mm},\end{equation}

from the stability conditions, ${\tau _c} = 9.18R{a^{ - 2/3}}$, ${k_c}( = 2{\rm \pi} d/{\lambda _c}) = 0.182R{a^{1/3}}$ and ${\tau _d} = 6{\tau _c}$. These scaling relationships are similar to those for the constant viscosity system, i.e. $\varGamma = 0$, given in (5.1a). Furthermore, these numerical values are independent of the initial liquid layer depth d, if $d \ge 1\;\textrm{cm}$. Even though we do not simulate the high Rayleigh number case due to the convergence issue of the present simulation, the scaling relations obtained in the present simulation can be used to explain the real experiment.

Table 3. Damköhler number, Rayleigh number, and the conversion factor for the dimensionless time to the real time based on the physical properties of NaCl–water system.

6. Conclusions

The formation of dissolution pattern on a horizontal salt body was investigated theoretically and numerically by considering the dissolution reaction and the diffusive and convective mass transfer at the dissolving interface. Such complex mass and momentum transport phenomena coupled with a dissolution reaction were modelled by introducing physical parameters and dimensionless numbers, in this case the Rayleigh number, $Ra$, the Damköhler number, $Da$, the saturated solute concentration, ${C_{sat}}$, the expansion factor, $\chi $, and the log-viscosity variation parameter, $\varGamma $. These were then solved using a linear stability analysis and a numerical simulation. As the dissolution of a solid salt proceeded, the level of the solution rose continuously and the solute became concentrated near the interface of the salt and the solution. As the level of the solution increased, the concentration gradient of solute resulted in an adverse density gradient. Due to the adverse density gradient, gravitational instability developed after a certain diffusion-dominant regime. Simultaneously, the instability motion enhanced the dissolution rate and pattern formation on the dissolving surface. The interplay of $Ra$, $Da$, ${C_{sat}}$, $\chi $ and $\varGamma $ on the onset of instability motion was predicted by a linear stability analysis, showing the same trend as the numerical simulation outcomes, suggesting scaling relationships for the pattern formation process. As a case study, pattern formation on a dissolving NaCl surface was numerically simulated, proving that the pattern formation in this case is mainly determined by the coupling of the dissolution rate and diffusive or convective mass transfer (the average height change, $R{a^{1/3}}\Delta {h_{avg}}\sim {(R{a^{2/3}}\tau )^{0.86}}$, for the convective mass transfer regime). The present scaling for the convection regime is slightly different from the experimentally observed one in the quasi-steady state regime, $R{a^{1/3}}\Delta {h_{avg}}\sim (R{a^{2/3}}\tau )$. Correspondingly, a 3-D numerical simulation visualized the initial stage of pattern formation on a dissolving surface.

This study elucidated the effects of $Ra$, ${C_{sat}}$, $\chi $ and $\varGamma $ on the nature of dissolution in the intermediate Da regime, i.e. ${10^{ - 2}} \le Da \le 100$, extending our understanding of the previous studies on dissolution in a horizontal geometry. Both the dissolution reaction and diffusive mass transfer determining Da play important roles in the dissolution-driven convective instability and pattern formation outcomes.