2.1. Reaction kinetics

We study a two-layer system consisting of miscible solutions initially separated by a horizontal contact plane (figure 1). The upper layer of the system is always of lower density than the lower one, which allows us to exclude the development of a global RT instability from the very beginning of the study. We assume that nitric acid $A$ (it can initially be either in the upper or lower layer) diffuses to react with sodium hydroxide $B$ to form their salt $C$ under the production of water. This process is accompanied by significant heat release $Q$. The standard enthalpy of the reaction is $57\ \textrm {kJ}\ \textrm {mol}^{-1}$. Such a neutralization reaction can be described by the simplified equation

(2.1)\begin{equation} {A} + {B} \to {C} + {Q}, \end{equation}

with the reaction rate characterized by the constant $K$. Thus, a second-order exothermic neutralization reaction defined by (2.1) has a comparatively simple, albeit nonlinear, kinetics.

Figure 1. Geometrical configuration of the two-layer miscible system and coordinate axes.

2.2. Hele-Shaw flow approximation

The Hele-Shaw cell, since its description by Hele Shaw (Reference Hele Shaw1898), has proven to be a useful instrument to study quasi-two-dimensional flows. These flows develop between two parallel flat plates separated by an infinitesimally small gap. As is well known, Hele Shaw's theory for single-phase flow results in the depth-averaged equations for pressure and two-dimensional velocity fields. In the simplest case, the Hele-Shaw approximation reduces the equation of motion to Darcy's law formulated initially for the fluid filtration through a porous medium. In this limit, the flow in the Hele-Show cell is completely analogous to the fluid filtration through a porous medium. The majority of works, which we cited above, have used the Darcy model to describe the phenomena in the system under consideration (Almarcha et al. Reference Almarcha, Trevelyan, Riolfo, Zalts, El Hasi, D'Onofrio and De Wit2010; Hejazi & Azaiez Reference Hejazi and Azaiez2012; Tsuji & Müller Reference Tsuji and Müller2012; Lemaigre et al. Reference Lemaigre, Budroni, Riolfo, Grosfils and De Wit2013; Kim Reference Kim2014; Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2015). It should be noted that the Darcy equation for the description of fluid flow in porous media is valid under the condition that the mean velocity varies slowly in space (Zeng, Yortsos & Salin Reference Zeng, Yortsos and Salin2003). In many cases, however, this condition is not met. For example, this may be due to a sharp change in fluid properties, as in the context of the RT problem Martin, Rakotomalala & Salin (Reference Martin, Rakotomalala and Salin2002), or in viscous fingering (Martin et al. Reference Martin, Rakotomalala, Talon and Salin2011), or due to the presence of the interface in the immiscible system, which requires the correct formulation of the boundary conditions for the velocity (Bratsun & De Wit Reference Bratsun and De Wit2004). The simplest way to correct this inadequacy in the model is to add the term suggested by Brinkman (Reference Brinkman1947). A more complicated correction implies taking into account the inertial term in a special form in the equation of motion, as was suggested by Ruyer-Quil (Reference Ruyer-Quil2001) and Martin et al. (Reference Martin, Rakotomalala and Salin2002).

The system we consider consists of two vertical solid plates placed parallel to each other with a small distance between them. Let $x$, $z$ be the directions parallel to the flat plates, and $y$ the perpendicular direction, with $h$ the gap between the plates at $y=\pm h/2$ (see figure 1). For the sake of clarity, suppose that the $z$-axis is anti-directed to gravity. The resulting closed cavity is defined by $0\leqslant x\leqslant L$ and $-H\leqslant z\leqslant H$.

We start with a three-dimensional Navier–Stokes equation, coupled with a continuity equation, implying that the fluid is incompressible

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

(2.3)\begin{gather}{\hat{\rho}}\left(\frac{\partial \boldsymbol{u}}{\partial t}+{\boldsymbol{u}}\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}\right) ={-} \boldsymbol{\nabla} p + \eta\nabla^2 {\boldsymbol{u}} + {\hat{\rho}}{\boldsymbol{g}} . \end{gather}

Here, ${\boldsymbol {u}}: (u_x, u_y, u_z)$ is the velocity, $\hat {\rho }$ is the solution density and $p$ is the pressure; $\eta$ and ${\boldsymbol {g}}$ stand for the dynamic viscosity and the acceleration vector due to gravity, respectively. Under the assumption that the fluid adheres to the solid walls, the boundary conditions for the velocity at the solid plates are

(2.4a,b)\begin{equation} y={\pm} \frac{h}{2} :\quad {\boldsymbol{u}} = 0. \end{equation}

We further assume that the gap width $h$ is small enough so that the fluid flow may be considered as quasi-two-dimensional, i.e. a Hele-Shaw approximation is applicable. Taking into account the boundary condition (2.4a,b), the velocity can be approximated by the following functions:

(2.5)\begin{equation} \left.\begin{gathered} u_x (x,y,z) = \frac{3}{2}\left(1-\frac{4y^2}{h^2}\right) U(x,z),\\ u_y (x,y,z) = 0,\\ u_z (x,y,z) = \frac{3}{2}\left(1-\frac{4y^2}{h^2}\right) V(x,z), \end{gathered}\right\} \end{equation}

where ${\boldsymbol {v}}: (U, V)$ is the two-component velocity field obtained after averaging the velocity $\boldsymbol {u}$ across the gap.

The approximations (2.5) should then be substituted into the Navier–Stokes equation (2.3) and averaged across the gap

(2.6)\begin{equation} \langle\cdots\rangle = \frac{1}{h}\int_{-{h}/{2}}^{{h}/{2}} \ldots\,\textrm{d}y . \end{equation}

As a result, we obtain the motion equation written in the Hele-Shaw approximation

(2.7)\begin{equation} \rho_0\left(\frac{\partial{\boldsymbol{v}}}{\partial t}+\frac{6}{5}{\boldsymbol{v}} \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{v}}\right) ={-} \boldsymbol{\nabla} p + \eta\Delta{\boldsymbol{v}} - \frac{12\eta}{h^2}{\boldsymbol{v}} +\langle\rho\rangle{\boldsymbol{g}}, \end{equation}

where $\rho _0$ is the density of the solvent and $\langle \rho \rangle$ is the average medium density, which will be determined below. In addition to the correction factor $6/5$ (Ruyer-Quil Reference Ruyer-Quil2001), (2.7) differs from the standard Navier–Stokes equation by the term proportional to the velocity. One can interpret this term as the average friction force due to the presence of the plates, and it is analogous to Darcy's law for the porous medium.

Let us estimate whether it is necessary to use the Hele-Shaw approximation in the form (2.7), which includes the Brinkman term and the inertial term, in this problem. We define the characteristic size of the chemoconvective structure as $H$, which can be estimated as $\sqrt {\tau _R D_{a0}}$, where $\tau _R=1/KA_0$ stands for the characteristic reaction time and $D_{a0}$ is the table value of the diffusion coefficient of acid. Then we scale the equation (2.7) choosing the following characteristic units: for time $\tau _R$, for length $H$, for velocity $H/\tau _R$, for pressure $12\rho _0 \nu D_{a0}/h^2$ and obtain

(2.8)\begin{equation} \frac{\varepsilon}{Sc}\left(\frac{\partial{\boldsymbol{v}}}{\partial t}+\frac{6}{5}{\boldsymbol{v}} \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{v}}\right) ={-} \boldsymbol{\nabla} p + \varepsilon\Delta{\boldsymbol{v}} - {\boldsymbol{v}} +\varepsilon Sc Ga\frac{\langle\rho\rangle}{\rho_0}{\boldsymbol{z}} , \end{equation}

where $Sc$ is the Schmidt number, $Ga$ is the Galileo number and

(2.9)\begin{equation} \varepsilon = \frac{h^2}{12H^2} = \frac{h^2\sqrt{KA_0}}{12\sqrt{D_{a0}}} \end{equation}

stands for the parameter determining the relative contribution of the Brinkman term to the motion equation.

In addition, we can evaluate the applicability of the Hele-Shaw approximation by reasoning in terms of the convective instability itself, regardless of the characteristic size of the cavity, which largely imposes its scale on physical processes. In our problem, there are two fundamental types of instabilities: RT convection and diffusive instability. Let us estimate the characteristic wavelength of each instability expressed in terms of diffusion, viscosity and buoyancy. An estimate of the characteristic length of diffusive types of instabilities can be obtained using the formula derived by Stern (Reference Stern1960). Stern considered doubly stratified water: warm and salty above cold and fresh, that is typical for the ocean stratification, and was the first to explain the mechanism of the onset of double-diffusion instability. Since, in our work, the role of heat is played by a rapidly diffusing acid, we must replace the characteristics of the thermal field in the Stern formula with the characteristics of the acid concentration field. Then, we obtain the estimate for an aqueous solution of nitric acid

(2.10)\begin{equation} L_{DD} = {\rm \pi}\left(\frac{4\nu D_{a0}}{g\beta_{a}|G_{a}|}\right)^{1/4} \approx 0.7\ \textrm{mm}, \end{equation}

where $\beta _a$ is the acid expansion coefficient and $G_a$ is the characteristic vertical gradient of acid concentration in the domain where the instability occurs. To obtain (2.10), we have used the data presented in table 1 of Part 1 (Mizev et al. Reference Mizev, Mosheva and Bratsun2021).

To estimate the characteristic size of the RT convection at the very beginning of the evolution, we can use the formula suggested by Martin et al. (Reference Martin, Rakotomalala and Salin2002). The authors have introduced the characteristic length $L_{RT}$ for the RT instability developing in miscible systems. For our case, we get the following estimate:

(2.11)\begin{equation} L_{RT} = \left(\frac{2\rho_0 \nu D_{a0}}{g \Delta\rho}\right)^{1/3} \approx 0.1\ \textrm{mm}, \end{equation}

where $\Delta \rho$ stands for the characteristic density difference giving rise to the RT instability.

We can learn from the estimates (2.10) and (2.11) that the characteristic instability length is either of the same order of magnitude as the Hele-Shaw cell gap ($h = 1\ \textrm {mm}$) used in the experiment or even smaller. This serves as another justification for using the Navier–Stokes–Darcy model equation (2.7) instead of the simple Darcy model. What is important is that the use of the averaged two-dimensional Navier–Stokes–Darcy equation (2.7) allows us to unify all limiting cases. The Darcy model for the Hele-Shaw cell description is valid when the cell gap is small compared to this characteristic reaction–diffusion length ($\varepsilon \ll 1$), whereas the case of fully three-dimensional flows governed by the Navier–Stokes equation is obtained in the opposite limit ($\varepsilon \gg 1$). The use of the equation (2.7) allows us to recover these two limits and to give a good approximation in the intermediate range of cell thicknesses (Ruyer-Quil Reference Ruyer-Quil2001; Martin et al. Reference Martin, Rakotomalala and Salin2002).

We can learn from the experimental observations (Mizev et al. Reference Mizev, Mosheva and Bratsun2021) that, at the very beginning of evolution (at $t<200\ \textrm {s}$ in dimensionless units), the characteristic size of the structures can be small enough and the diffusion represented by the Brinkman term in (2.7) can play an important role in the system dynamics. At later stages of evolution ($t>200\ \textrm {s}$), Darcy's law can be used to model flows. In fact, the situation is ambiguous, since chemoconvective cells form for a rather long time. In the case of the SW mode, the processes are much faster, but the characteristic size of convective structures is larger. Summarizing all of this, we can state that the use of the full three-dimensional Navier–Stokes equations in this problem is redundant, but the use of the Brinkman term is desirable.

All other physical fields should be averaged in the same style as the velocity.

2.3. Boussinesq approximation

The Boussinesq approximation is commonly applied in convection theory to describe buoyancy-driven flows. This approximation embodies two essential ideas. First, any fluctuations in density that occur with the onset of fluid motion are produced mainly by either thermal effects (in the case of thermally induced buoyancy) or by concentration effects (in the case of solutally induced buoyancy) rather than by pressure effects. The second idea is that all variations in fluid properties, except density, may be neglected. Furthermore, this approximation ignores density differences except where they appear in the terms responsible for the convective buoyancy force. The Boussinesq approach is justified for the problems in which ‘weak’ convection exists in the cavity on the laboratory scale, and the density variations caused by thermal or concentration expansion are relatively small.

In fact, we have already implicitly used the Boussinesq approach when we replaced the density function $\rho$ on the left side of the motion equation (2.7) with a constant value $\rho _0$ of the solvent density. To complete the work, one needs to expand the density $\langle \rho \rangle$ on the right-hand side of (2.7) as a power series of temperature and concentrations, retaining only linear terms

(2.12)\begin{equation} \langle\rho\rangle = \rho_0 + \rho = \rho_0(1 -\beta T +\beta_a A +\beta_b B +\beta_c C) , \end{equation}

where $\beta$ and $\beta _i$ are the thermal and solutal expansion coefficients, respectively. In (2.12), we take into account the fact that all the dissolved substances are heavier than water. Let us evaluate the contribution of thermal and solutal effects to the buoyancy force by composing the ratio of the first two terms of the first order in the expansion (2.12)

(2.13)\begin{equation} \frac{\beta\varTheta^*}{\beta_A A^{*}} = Le\frac{\beta|Q|D_{a0}}{\beta_a \kappa}\approx 0.08 , \end{equation}

where $D_{a0}$ is the table value of the diffusion coefficient of acid, $\kappa$ is the thermal conductivity of solvent and $\varTheta ^{*}$ and $A^{*}$ stand for the characteristic values of the temperature and acid concentration, respectively; $Le$ stands for the Lewis number defined as the ratio of thermal diffusivity to mass diffusivity. To get an estimate (2.13), we have used data for water and nitric acid: $Le=42$, $\beta = 0.207\times 10^{-3}\ \textrm {K}^{-1}$, $\beta _a$ = $3.34\times 10^{-2}\ \textrm {l}\ \textrm {mol}^{-1}$, $D_{a0} = 3.15\times 10^{-5}\ \textrm {cm}^2\ \textrm {s}^{-1}$, $\kappa = 0.6\ \textrm {J}\ \textrm {s}^{-1}\ \textrm {m}^{-1}\ \textrm {K}^{-1}$.

In addition to the evaluation (2.13), one can notice that the thick walls of the experimental Hele-Shaw cell are made from glass, which transmits significant heat because the thermal conductivity coefficients of water and glass are nearly the same. In comparison with the concentration effects, the thermal effects can be controlled to a greater extent during the experiment. Thermally insulated walls enhance the role of heat, while perfectly conductive walls make the heat effect negligible.

Finally, one can learn from (2.13) that the thermally induced buoyancy is an order of magnitude weaker than the buoyancy force caused by the solute effect. For this reason, in what follows, we consider the problem as isothermal.

2.4. Governing equations

As follows from the experimental observations presented in Part 1 of this work (Mizev et al. Reference Mizev, Mosheva and Bratsun2021), the processes that occur during the neutralization of a strong acid with alkali metal hydroxides are very similar. Therefore, we present the theoretical study for one characteristic pair of reactants: nitric acid and sodium hydroxide. So, let two miscible fluids fill the Hele-Shaw cell with a gap depth $h$, and let them be spatially separated at the very beginning by a contact surface $z=0$. The upper (or lower) layer is an aqueous solution of $\textrm {HNO}_3$, and the lower (or upper) layer is an aqueous solution of $\textrm {NaOH}$ (figure 2). The initial values of reactant concentrations are $A_0$ and $B_0$. These two parameters are key to the problem since they are the only ones we can easily tune in the experiments with a fixed pair of reactants. In all cases considered below, the initial concentrations of reactants are selected so that the stratification of the system at the very beginning is statically stable. Thus, we exclude the occurrence of the RT instability.

Figure 2. Instantaneous density profiles $\rho ^0(z,t)$ calculated for successive times in the case of pure diffusion (a), reaction–diffusion without the CDD effect (b) and reaction–diffusion influenced by the CDD effect ($c$). The position of the initial contact surface of the solutions is determined by the line $z = 0$. The initial concentrations are $\gamma _a=0.667$, $\gamma _b=0.7$.

We use $h$ as the measurement unit for length and $h^2/D_{a0}$ for time, $D_{a0}/h$ for velocity, $\rho _0\nu D_{a0}/h^2$ for pressure and $A_{lim}$ for concentration, where $A_{lim}$ is the maximum concentration at which the linear law of CDD is observed for nitric acid. The incompressibility constraint (2.2) can be satisfied automatically by defining a streamfunction $\varPsi$ such that

(2.14a,b)\begin{equation} U ={-} \frac{\partial\varPsi}{\partial z},\quad V = \frac{\partial\varPsi}{\partial x} . \end{equation}

In what follows, we use the vorticity–streamfunction formulation of the governing equations. Thus, the set of reaction–diffusion–advection equations written in the Boussinesq and Hele-Shaw approximations is as follows:

(2.15)\begin{gather} \varPhi ={-} \nabla^2 \varPsi , \end{gather}

(2.16)\begin{gather}\frac{1}{Sc}\left(\frac{\partial\varPhi}{\partial t}+ \frac{6}{5}\frac{\partial\varPsi}{\partial z}\frac{\partial\varPhi}{\partial x} - \frac{6}{5}\frac{\partial\varPsi}{\partial x}\frac{\partial\varPhi}{\partial z}\right) = \nabla^2 \varPhi - 12\varPhi - R_{a}\frac{\partial A}{\partial x} - R_{b}\frac{\partial B}{\partial x} - R_{c}\frac{\partial C}{\partial x}, \end{gather}

(2.17)\begin{gather}\frac{\partial A}{\partial t} + \frac{\partial\varPsi}{\partial z} \frac{\partial A}{\partial x} - \frac{\partial\varPsi}{\partial x} \frac{\partial A}{\partial z} = \boldsymbol{\nabla} (D_a\boldsymbol{\nabla} A) - Da A B, \end{gather}

(2.18)\begin{gather}\frac{\partial B}{\partial t} + \frac{\partial\varPsi}{\partial z} \frac{\partial B}{\partial x} - \frac{\partial\varPsi}{\partial x} \frac{\partial B}{\partial z} = \boldsymbol{\nabla} (D_b\boldsymbol{\nabla} B) - Da A B , \end{gather}

(2.19)\begin{gather}\frac{\partial C}{\partial t} + \frac{\partial\varPsi}{\partial z} \frac{\partial C}{\partial x} - \frac{\partial\varPsi}{\partial x} \frac{\partial C}{\partial z} = \boldsymbol{\nabla} (D_c\boldsymbol{\nabla} C) + Da A B , \end{gather}

where (2.15) is the definition of the vorticity $\varPhi$. The diffusion terms in (2.17)–(2.19) have been written in the most general form taking into account the CDD effect (Crank Reference Crank1975).

Equations (2.15)–(2.19) should be supplemented by the boundary conditions

(2.20a,b)\begin{gather} x = 0,L:\quad\varPsi = 0, \quad\frac{\partial\varPsi}{\partial x} = 0 , \quad\frac{\partial A}{\partial x}=0 , \quad\frac{\partial B}{\partial x}=0 , \quad\frac{\partial C}{\partial x}=0 , \end{gather}

(2.21a,b)\begin{gather}z ={\pm} L:\quad\varPsi = 0, \quad\frac{\partial\varPsi}{\partial z} = 0 , \quad\frac{\partial A}{\partial z}=0 , \quad\frac{\partial B}{\partial z}=0 , \quad\frac{\partial C}{\partial z}=0 . \end{gather}

If the acid solution is on the top, the initial conditions at $t=0$ are written as

(2.22a,b)\begin{gather} z \leqslant 0: \quad \varPsi = 0,\quad A = 0, \quad B=\gamma_b , \end{gather}

(2.23a,b)\begin{gather}z > 0:\quad\varPsi = 0,\quad A = \gamma_a, \quad B=0 , \end{gather}

where $\gamma _a$ and $\gamma _b$ stand for

(2.24a,b)\begin{equation} \gamma_a = \frac{A_0}{A_{lim}} , \quad \gamma_b = \frac{B_0}{A_{lim}} , \end{equation}

respectively.

If the base solution is on the top, the initial conditions are

(2.25a,b)\begin{gather} z \leqslant 0: \quad \varPsi =0,\quad A = \gamma_a, \quad B=0, \end{gather}

(2.26a,b)\begin{gather}z > 0:\quad\varPsi = 0,\quad A = 0, \quad B=\gamma_b. \end{gather}

The set of equations (2.15)–(2.19) has two dimensionless parameters

(2.27a,b)\begin{equation} Sc = \frac{\nu}{D_{a0}}, \quad Da = \frac{K A_{lim}h^2}{D_{a0}} , \end{equation}

the Schmidt number and the Damköhler number, respectively. These two numbers are the ratio of either the viscous diffusion rate or reaction rate to the mass diffusion rate, respectively. The evaluation of the Schmidt number for nitric acid gives the value $Sc = 317$. As for the Damköhler number, in this paper, we accept a comparatively high value, which makes the reaction frontal: $Da = 10^3$. Note that the above parameters (2.27a,b) are determined by the reaction kinetics, and they cannot be tuned in the experiments with a particular pair of reactants.

In addition, the set of Rayleigh numbers appears in (2.16)

(2.28a–c)\begin{equation} R_a = \frac{g\beta_a A_{lim}h^3 }{\nu D_{a0}} , \quad R_b = \frac{g\beta_b A_{lim}h^3 }{\nu D_{a0}} , \quad R_c = \frac{g\beta_c A_{lim}h^3 }{\nu D_{a0}} , \end{equation}

characterizing the contribution of each water dissolved substance to density variations. The values of the parameters (2.28a–c) can also be estimated with the table values: $R_{a}=3.18\times 10^5$, $R_{b}=3.82\times 10^5$, $R_{c}=5.1\times 10^5$.

In what follows, it is convenient to define a variable based on the expansion (2.12), which makes sense of the addition to the density of a pure solvent due to dissolved substances

(2.29)\begin{equation} {\rho} (x,z,t) = A(x,z,t) + \frac{R_b}{R_a} B(x,z,t) + \frac{R_c}{R_a} C(x,z,t) . \end{equation}

Finally, the only parameters which can be manipulated by an experimenter are the initial values of the reactant concentrations $\gamma _a$ and $\gamma _b$.

Thus, by obtaining a full set of equations (2.15)–(2.19) with the boundary conditions (2.20a,b), (2.21a,b) and the initial conditions in the form of either (2.22a,b), (2.23a,b) or (2.25a,b), (2.26a,b), we have formulated the problem. Note that the reaction-induced buoyancy number $K_{\rho }$ does not automatically appear in the governing equations (2.15)–(2.19) after they have been converted to a dimensionless form. We will define this parameter below.

2.5. Concentration-dependent diffusion

The CDD effect is key to this work. To evaluate the diffusion formulas for the pair $\textrm {HNO}_3/\textrm {NaOH}$, we have brought together all the known experimental data (Bratsun et al. Reference Bratsun, Kostarev, Mizev and Mosheva2015). Note that the data on the concentration dependence of the diffusion coefficients have appeared to be fragmentary and incomplete for most substances. Therefore, we were forced to conduct our experiments to measure the diffusion coefficients (see data presented in Mizev et al. Reference Mizev, Mosheva and Bratsun2021). We assume for simplicity that the data set falls in the experimentally interesting range of concentration on a straight line $f(X) = a + bX$, where $X$ stands for concentration, $a$ and $b$ are some constants. We suppose here that each diffusion coefficient depends only on the concentration of its substance. These linear laws in a dimensionless form, then, are

(2.30)\begin{equation} \left.\begin{gathered} D_a (A) \approx 0.881 + 0.158 A , \\ D_b (B) \approx 0.594 - 0.087 B ,\\ D_c (C) \approx 0.478 - 0.284 C . \end{gathered}\right\} \end{equation}

Equations (2.30) perfectly approximate the diffusive properties within the concentration range from $0.1$ up to $A_{lim} = 3\ \textrm {mol}\ \textrm {l}^{-1}$. The dimensionless table values for diffusion coefficients are as follows:

(2.31a–c)\begin{equation} D_a^{tab}=1.0,\quad D_b^{tab}=0.68,\quad D_c^{tab}=0.5 , \end{equation}

where the coefficients were scaled using $D_{a0}$.

One can see that the results of (2.30) for the limiting case of zero concentration do not coincide with the table values (2.31a–c). This is explained by the fact that the experimental dependencies of the diffusion coefficients at small concentrations are strongly nonlinear. Since the neutralization reaction under ideal mixing runs fast, we can neglect these nonlinear effects at low concentrations.

3.1. Potential barrier and well

The problem (2.15)–(2.26a,b), (2.30) has an unsteady solution, which describes the reaction–diffusion processes. The fluid remains at mechanical equilibrium. We will consider this solution as a base-state solution. We assume that the fluid velocity equals zero in (2.15)–(2.19) and that the concentration fields depend only on the vertical coordinate and time: $A^0(z,t)$, $B^0(z,t)$, $C^0(z,t)$. The resulting time-dependent nonlinear equations

(3.1)\begin{gather} \frac{\partial A^0}{\partial t} = D_{a}(A^0)\frac{\partial^2 A^0}{\partial z^2} + \frac{\textrm{d} D_{a}(A^0)}{\textrm{d}A^0}\left(\frac{\partial A^0}{\partial z}\right)^2 - Da A^0 B^0 , \end{gather}

(3.2)\begin{gather}\frac{\partial B^0}{\partial t} = D_{b}(B^0)\frac{\partial^2 B^0}{\partial z^2} + \frac{\textrm{d} D_{b}(B^0)}{\textrm{d}B^0}\left(\frac{\partial B^0}{\partial z}\right)^2 - Da A^0 B^0 , \end{gather}

(3.3)\begin{gather}\frac{\partial C^0}{\partial t} = D_{c}(C^0)\frac{\partial^2 C^0}{\partial z^2} + \frac{\textrm{d} D_{c}(C^0)}{\textrm{d}C^0}\left(\frac{\partial C^0}{\partial z}\right)^2 + Da A^0 B^0 \end{gather}

should be complemented with the boundary and initial conditions (2.20a,b)–(2.26a,b) and the formulas for CDD (2.30). This problem has no analytical solution and can only be solved numerically.

Figure 2 presents the base-state profiles of the density $\rho ^0 (z,t)$ for four consecutive times $t=0,2,5,10$. We consider three fundamentally different cases. Figure 2(a) illustrates the situation when, the solutions being brought into contact, do not react ($Da=0$). In this case, one can observe a typical profile of the DLC instability, which occurs when the system is initially statically stable, and the fastest diffusing solute is in the upper layer. One can see here that a sharp drop in the density of the solutions given at the very beginning is gradually smoothed out since diffusion is a typical relaxation process.

Figure 2(b) shows the dynamics of the density profile starting from the same initial conditions, but now the solutions of $\textrm {HNO}_3$ and $\textrm {NaOH}$ react according to the equation (2.1). The CDD effect (2.30) is here not taken into account. Instead, the equations are integrated with constant tabular values (2.31a–c) for the diffusion coefficients. One can see that the situation does not change qualitatively, but the profiles become non-symmetric concerning the $z$-axis. However, the main feature of the density profile, to spread out over time under the influence of the relaxation process, has been preserved.

Finally, figure 2(c) shows the density evolution taking into account the CDD effect. A qualitative change here is the appearance of a potential barrier on the density profile. For this reason, the density profile has two minima above and below the reaction. It is interesting to note that the height and width of the potential barrier vary very little with time. This implies that the barrier is the result of an exact balance between CDD and the nonlinear chemical reaction.

Figure 3 shows the relative contributions of acid, base and salt to the total density $\rho ^0(z,t)$ of the medium. One can see that the reaction product makes a decisive contribution to the emergence of a potential barrier. Due to the barrier, which successfully resists the relaxation process, a potential well that is statically quasi-stable in time also arises in the system. Since gravity is anti-directed to the $z$-axis, only the potential well can exist (the region to the left of the barrier in figure 3). The bottom of the potential well is marked with an open square. In this case, the potential barrier prevents the low-density area from floating up under the action of the Archimedes force. The low-density region to the right of the barrier is statically unstable and can float under the action of a buoyancy force. It is interesting to note that the well's depth almost does not change with time, yet its width increases time as the left wall is gradually smoothed out due to diffusion (see figure 2c).

Figure 3. Base-state profiles of acid $A^0$, base $B^0$ and their salt $C^0$ showing contributions of species to the total density $\rho ^0$. The initial concentrations are $\gamma _a=0.667$, $\gamma _b=0.55$. All profiles are shown for time $t=5$. The position of the potential barrier is indicated by the circle.

Thus, the neutralization reaction, coupled with the CDD, can result in a potential well and maintain it for a long time in a quasi-steady state. As we will show below, this ability provides completely new opportunities for pattern formation scenarios in the system.

3.2. Stability map based on the base-state density profiles

We found that the stability map of the principal reaction–diffusion–convection modes observed experimentally (Mizev et al. Reference Mizev, Mosheva and Bratsun2021) can be constructed by analysing changes in the base-state density profile $\rho ^0 (z,t)$ calculated using (3.1)–(3.3). Figure 4 presents this map in the parameter plane of the initial values for the concentrations of acid $\gamma _a$ and base $\gamma _b$ defined by (2.24a,b).

Figure 4. Stability map constructed by inspection of variations of the base-state density profile in the ($\gamma _a, \gamma _b$) parameter space. Abbreviations DLC, SW and CDD denote the diffusive layer convection ($1$), shock wave ($3$) and convection of CDD ($2$), respectively. The characteristic cross-section $\gamma _a =0.667$ of the stability map is indicated by a red straight line and discussed in figures 5 and 6 in more detail. Experimental interferograms (Mizev et al. Reference Mizev, Mosheva and Bratsun2021) are shown alongside to illustrate the main modes indicated in the map.

The relation (2.29) implies that the weight of an elementary volume of liquid depends both on the structure of the solute's molecule and on the amount of solute. Therefore, the condition

(3.4)\begin{equation} \gamma_b = \frac{\beta_a}{\beta_b}\gamma_a = \frac{R_a}{R_b}\gamma_a \end{equation}

on the stability map gives a line of equal density for the upper and lower layers, that is, an isopycnal line (thick solid line in figure 4). If the parameters are taken above the isopycnal line, then the base solution is heavier than the acid solution. Since we do not consider here the RT instability, this means that a less dense fluid lies above a denser fluid in a gravity field. Therefore, the intersection of the isopycnal line (3.4) means that the solutions change places.

After bringing the solutions of $\textrm {HNO}_3$ and $\textrm {NaOH}$ into contact, the reaction–diffusion processes transform the density profile and may cause potentially unstable conditions for the system under the gravity force. Density transformations under different initial concentrations are shown in figure 5. This figure presents typical density profiles that correspond to a vertical slice indicated in figure 4, which includes the main areas above and below the isopycnal line. Let us discuss the characteristic areas of the stability map in more detail. The upper left part of the map corresponds to the parameters at which the DLC instability arises (the solution with a faster diffusing substance is at the top). As a rule, with no reaction, it is a formation of two areas with weak convective movements divided by a diffusion zone where the fluid remains motionless. The zones above and below are characterized by the system of fingers which propagate symmetrically up and down from the mixing zone (Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011). It is known that the running reaction could violate the symmetry of the fingering process, giving preference to fingers that extend up or down (Almarcha et al. Reference Almarcha, Trevelyan, Riolfo, Zalts, El Hasi, D'Onofrio and De Wit2010). However, the effect of the reaction on the mass transfer processes in the system can be much more radical (Bratsun et al. Reference Bratsun, Kostarev, Mizev and Mosheva2015, Reference Bratsun, Mizev, Mosheva and Kostarev2017).

Figure 5. Instantaneous density profiles $\rho ^0(z,t)$ for the reaction–diffusion base state calculated for different values of the control parameters $\gamma _a$ and $\gamma _b$ marked on the stability map (see figure 4). Red horizontal lines are drawn for the convenience of comparing the reaction zone and upper layer densities. The position of the initial contact surface of the solutions is determined by $z = 0$. All profiles are shown for time $t=5$.

Moving down along the red line $\gamma _a=0.667$ shown in figure 4, we first cross the point $a$, which corresponds to the appearance of a new local maximum on the density curve (figure 5a). The reason is that the diffusion rate of salt decreases with the concentration increase (see (2.30)), which results in local accumulation of the reaction product near the initial contact surface of the solutions. Also, since the salt is quite heavy, its contribution to the density (2.29) is significant. Thus, the density maximum indicates where the reaction product accumulates in the system and serves as an indicator for the reaction front location. The appearance of a potential well, or in other words ‘a density pocket’, in figure 5(a) does not mean the immediate development of convection there. For this, the local Rayleigh number must exceed its critical value. Let us define this parameter as follows:

(3.5)\begin{equation} R_{local} = \frac{1}{R_{local}^{cr}}\frac{g\beta_c N d^4}{16 \nu D_{c0}}, \end{equation}

where $D_{c0}$ is the table value of the diffusion coefficient of salt, $N$ is the characteristic linear gradient of salt concentration in a potential well and $d$ is the width of the potential density well. The semi-width $d/2$ is chosen as the characteristic size because the instability condition is fulfilled only in one half of the potential well. It can be calculated as

(3.6)\begin{equation} \frac{d}{2}=|z_{bar}-z_{cdd}| \end{equation}

(see figure 3 for details). The parameter $R_{local}$ is normalized by its critical value $R_{local}^{cr}$, which could be either calculated within a linear stability analysis or taken from a textbook on the solutal Rayleigh–Bénard problem (Gershuni & Zhukhovitskii Reference Gershuni and Zhukhovitskii1976). In the definition (3.5), we also take into account the fact that the main contribution to the appearance of the potential well is made by the reaction product.

The variation of $R_{local}$ as a function of $\gamma _b$ and fixed $\gamma _a = 0.667$ is shown in figure 6(a). It is seen that the fluid loses its stability with respect to local perturbations limited to the density pocket at approximately $\gamma _b = 0.73$ (point $b$ in figure 4). The corresponding density profile is presented in figure 5(b). We have called this type of instability CDD convection, emphasizing the fact that, without the effect of CDD, convection in the potential well would not have been possible. To demonstrate this, figure 6(a) also presents the curve calculated without taking into account the CDD effect. Although the local Rayleigh number slightly exceeds the critical value of 1 near the isopycnal line itself, this instability still does not occur because, in this parameter range, the system becomes unstable to another type of perturbation (see below). In fact, the potential well, where the CDD instability occurs, is destroyed by a global bifurcation and begins to float upward together with the entire reaction zone.

Figure 6. Variation of the local Rayleigh number calculated inside a potential well (a) and the inverse of the reaction-induced buoyancy number $K_\rho$ (b) as a function of the dimensionless initial concentration of base $\gamma _b$. The corresponding cross-section of the stability map is given in figure 4. The critical value of the parameter in each case is indicated by a red dashed line (instability above the line).

By inspecting the critical density profile shown in figure 5(b), we can conclude the reaction zone density is still less than the density of the lower layer, but higher than that of the upper layer. This allows the DLC and CDD instabilities to coexist together while being separated by a thin diffusive layer.

Moving further down along the red line in figure 4, we cross the point $c$ ($\gamma _b\approx 0.58$), where the densities of the reaction zone and the upper layer become equal. In what follows, for convenience, we define the reaction-induced buoyancy number $K_{\rho }$ introduced in Mizev et al. (Reference Mizev, Mosheva and Bratsun2021) as

(3.7)\begin{align} K_{\rho} = \begin{cases} \dfrac{1 + A^0(z_{bar},t_c) + \dfrac{R_b}{R_a}B^0(z_{bar},t_c) + \dfrac{R_c}{R_a} C^0(z_{bar},t_c)}{1 + A^0(H,t_c)}, & \text{if acid is above base};\\ \dfrac{1 + A^0(z_{bar},t_c) + \dfrac{R_b}{R_a}B^0(z_{bar},t_c) + \dfrac{R_c}{R_a} C^0(z_{bar},t_c)}{1 + \dfrac{R_b}{R_a} B^0(H,t_c)}, & \text{if base is above acid}, \end{cases} \end{align}

where $z_{bar}$ stands for the position of the local maximum density in the reaction zone (it coincides with the position of the potential barrier, see figure 3), $t_c$ is the point in time when the density measurement occurs. Generally, the measurement should be carried out in the limit of the large time asymptotic solution $t_c\to \infty$. In practice, the parameter (3.7) was calculated for a sufficiently long time $t_c=10$.

In figure 6(b), we show the variation of $1/K_{\rho }$ as a function of $\gamma _b$ and fixed $\gamma _a = 0.667$. The density profile in figure 5(c) indicates that this is the point for global bifurcation in the system, at which an intense convective motion in the form of a shock wave occurs (Bratsun et al. Reference Bratsun, Mizev, Mosheva and Kostarev2017). It can be noted that the area limited by the curve $K_{\rho } = 1$ is a narrow band pressed to the isopycnal line. The area becomes narrower and fades away at both small and comparatively large values of the initial concentrations. The estimates show that at $\gamma _a > 1$ the effect disappears.

Moving further along the parameter plane, we pass over the isopycnal line at $\gamma _b=0.55$ (the point $d$ in figure 4) and enter the parameter area where the acid and base solutions change places. Figure 6(b) demonstrates that, below the isopycnal line, there is also a region of intense convective motion. The exit from this mode occurs at the point $e$ ($\gamma _b\approx 0.52$). The corresponding density profile is shown in figure 5(e).

It is important to show why the CDD convection does not occur if the acid and base change places in the system (below the isopycnal line in figure 4). Figure 5 clearly demonstrates that the bifurcation points for the onset of the CDD and SW instabilities practically coincide. Due to this, the CDD convection, which is principally local, does not have a chance to develop. Ultimately, this happens because of the asymmetry of the local maximum with respect to the up–down reflection, which is the result of different diffusion rates of $\textrm {HNO}_3$ and $\textrm {NaOH}$. The slope of the potential wall located on the side of the base solution is always greater since the diffusion rate of nitric acid in the water is higher (for example, compare the profiles shown in the figure 5c,e). Thus, chemoconvective cells cannot arise, although the potential well supported by the CDD effect and nonlinear reaction continues to exist between the points $e$ and $f$ in figure 4.

4.1. Analytical QSSA solution for the simplified problem

In our problem, an initially stable density profile evolves to a state with a potential well. One can see in figure 2(c) that the density profile changes rapidly at the very beginning and then slowly evolves, but does not change qualitatively. We can use this property employing a linear theory analytically within the QSSA approach.

Let us consider the following simplified system of equations:

(4.1)\begin{gather} \varPhi + \nabla^2 \varPsi = 0 , \end{gather}

(4.2)\begin{gather}\nabla^2 \varPhi - 12\varPhi - \frac{\partial \rho}{\partial x} = 0 , \end{gather}

(4.3)\begin{gather}\frac{\partial \rho}{\partial t} + \frac{\partial\varPsi}{\partial z}\frac{\partial\rho}{\partial x} - \frac{\partial\varPsi}{\partial x}\frac{\partial\rho}{\partial z} = \nabla^2 \rho , \end{gather}

where $\rho$ is the density of the medium. The variables in (4.1)–(4.3) were rescaled in such a way that the diffusion coefficient and Rayleigh number were removed from the equations (Hejazi et al. Reference Hejazi, Trevelyan, Azaiez and De Wit2010; Gandhi & Trevelyan Reference Gandhi and Trevelyan2014). We neglect the processes of reaction and nonlinear diffusion in (4.3). To simplify the analysis, we replace the reaction–diffusion equations generating a non-monotonic density profile by the initial state in the form of a rectangular potential well, which slowly (in terms of the development of hydrodynamic disturbances) expands with time

(4.4)\begin{equation} \rho^0 (\tau,z) = \left\{\begin{array}{@{}ll} 0 , & z > d(\tau)\\ -D , & 0\leqslant z \leqslant d(\tau) \\ 0 , & z<0\\ \end{array}\right. , \end{equation}

where $d$ stands for the width of the potential well changing over slow time $\tau$, $D$ is the constant well depth. The density profile constructed from a linear combination of step functions as in (4.4) allows us to obtain the closed-form analytical solution (Tan & Homsy Reference Tan and Homsy1986; Gandhi & Trevelyan Reference Gandhi and Trevelyan2014).

We decompose the small monotonic perturbations in the following way:

(4.5)\begin{equation} \begin{pmatrix} \varPhi(t,x,z)\\ \varPsi(t,x,z)\\ \rho (t,x,z) \end{pmatrix} = \begin{pmatrix} 0\\ 0\\ \rho^0 (\tau,z) \end{pmatrix} + \begin{pmatrix} {\varphi} (z)\\ {\psi} (z)\\ a (z) \end{pmatrix}\exp({\sigma t+ Ikx}), \end{equation}

where ${\varphi }, {\psi }, a$ are, respectively, the amplitudes of normal mode perturbations for the vorticity, streamfunction and density; $k$ is the wavenumber and $\sigma$ is the growth rate.

Substituting (4.5) into (4.1)–(4.3) and linearizing these equations near the base state (4.4), we obtain the following system of amplitude equations for the determination of critical perturbations:

(4.6)\begin{gather} \varphi + \frac{\textrm{d}^2\psi}{\textrm{d} z^2} - k^2\psi = 0 , \end{gather}

(4.7)\begin{gather}\frac{\textrm{d}^2\varphi}{\textrm{d} z^2} - (k^2+12)\varphi - k^2 a = 0 , \end{gather}

(4.8)\begin{gather}\sigma a = \frac{\textrm{d}^2 a}{\textrm{d} z^2} - k^2 a- \psi\frac{\textrm{d}\rho^0}{\textrm{d} z} . \end{gather}

Equations (4.6)–(4.8) should be supplemented by the condition that the solutions are continuous at the jumps

(4.9)\begin{equation} \left.\begin{gathered} z = 0, d(\tau):\quad \left[\psi\right]^{+}_{-} = 0,\quad \left[\varphi\right]^{+}_{-} = 0 ,\quad \left[a\right]^{+}_{-} = 0 ,\\ \left[\frac{\textrm{d}\psi}{\textrm{d}z}\right]^{+}_{-} = 0,\quad \left[\frac{\textrm{d}\varphi}{\textrm{d}z}\right]^{+}_{-} = 0,\quad \left[\frac{\textrm{d} a}{\textrm{d}z} + \psi\rho^0\right]^{+}_{-} = 0 . \end{gathered}\right\} \end{equation}

The sixth-order system of ordinary differential equations (4.6)–(4.8) can be analytically solved for three spatial ranges (inside and outside the potential well). The obtained solutions then are substituted into the jump conditions (4.9), yielding a linear system of twelve algebraic equations for the unknowns. We should require that the determinant $12\times 12$ coefficient matrix be equal to zero, which results in the dispersion equation

(4.10)\begin{align} \frac{24\sigma^2(\sigma-12)^2}{D^2k^2} &= \frac{k}{l}(\sigma-12)\left(1-\textrm{e}^{-(k+l)d}\right) - \frac{k\sigma}{12m}(\sigma-12)\left(1-\textrm{e}^{-(k+m)d}\right) \nonumber\\ &\quad -\frac{k^2\sigma}{ml}\left(1-\textrm{e}^{-(l+m)d}\right) + \frac{6 k^2}{l^2}\left(1-\textrm{e}^{{-}2ld}\right) +\frac{k^2\sigma^2}{24m^2}\left(1-\textrm{e}^{{-}2md}\right) \nonumber\\ &\quad +\frac{(\sigma-12)^2}{24}\left(1-\textrm{e}^{{-}2kd}\right), \end{align}

where $l^2\equiv k^2+\sigma$, $m^2\equiv k^2 + 12$.

The equation for the cutoff wavenumber $k$ is then obtained in the limit $\sigma \to 0$

(4.11)\begin{equation} d^2D^2 = \frac{2304 k^2 d^2}{F_1 +F_2\left(1-\textrm{e}^{-(k+m)d}\right)-F_3\,\textrm{e}^{{-}2kd}+F_4\left(1-\textrm{e}^{{-}2md}\right)}, \end{equation}

where

(4.12)\begin{equation} \left.\begin{gathered} F_1 = 1-\frac{1}{3}k^2 - \frac{12d-m}{36m}k^4,\quad F_2 = \frac{k^3}{3m}\left(1+dk - \frac{1}{6}k^2\right)\\ F_3 = 1+2dk +\frac{2d^2-1}{3}k^2 - \frac{H}{3}k^3+\frac{1}{36}k^4,\quad F_4 = \frac{k^6}{36(k^2+12)} . \end{gathered}\right\} \end{equation}

Let us determine the condition for the instability onset in the LW limit expanding the cutoff condition (4.11) for small $k$

(4.13)\begin{equation} \left({-}6d^2 +\frac{13824}{D^2}\right) k^2 + \left(2d - \frac{2}{m} + \frac{2}{m e^{md}} +4d^3\right) k^3 + \cdots.\end{equation}

To leading order in $k$, we obtain the important relation for the threshold

(4.14)\begin{equation} d(\tau)>\frac{48}{D} . \end{equation}

We can compare our result with the analytical result obtained by Tan & Homsy (Reference Tan and Homsy1986) and Gandhi & Trevelyan (Reference Gandhi and Trevelyan2014) for initial time $t=0$ in the case of the Darcy motion equation

(4.15)\begin{equation} d>\frac{4}{D} . \end{equation}

Our analysis (4.14) is applicable at any time. At the very beginning of evolution, the width of the potential well is extremely small; therefore, LW disturbances will not be excited. We should expect that the most dangerous disturbances will be short waves (SW). Comparing (4.14) and (4.15), we can conclude that taking into account the Brinkman term results in the fact that LW disturbances are successfully damped.

4.2. Numerical IVP solution

In our case, the instability always develops after some critical time and the IVP usage is more reasonable. Let us decompose the small monotonic perturbations in the following way:

(4.16)\begin{equation} \begin{pmatrix} \varPhi(t,x,z)\\ \varPsi(t,x,z)\\ A (t,x,z)\\ B (t,x,z)\\ C (t,x,z) \end{pmatrix} = \begin{pmatrix} 0\\ 0\\ A^0 (t,z)\\ B^0 (t,z)\\ C^0 (t,z) \end{pmatrix} + \begin{pmatrix} {\hat{\varphi}} (t,z)\\ {\hat{\psi}} (t,z)\\ a (t,z)\\ b (t,z)\\ c (t,z) \end{pmatrix} \textrm{e}^{Ikx}, \end{equation}

where ${\hat {\varphi }}, {\hat {\psi }}, a, b, c$ are, respectively, the amplitudes of normal mode perturbations for the vorticity, streamfunction, acid, base and salt concentrations while $k$ is their wavenumber.

Substituting (4.16) into (2.15)–(2.19) and linearizing these equations near the base state (3.1)–(3.3), we obtain the following system of time-dependent amplitude equations for the determination of critical perturbations:

(4.17)\begin{gather} \varphi + \frac{\partial^2\psi}{\partial z^2} - k^2\psi = 0,\end{gather}

(4.18)\begin{gather} \frac{1}{Sc}\frac{\partial\varphi}{\partial t} = \frac{\partial^2\varphi}{\partial z^2} - (k^2+12)\varphi - k^2 (Ra_a a + Ra_b b + Ra_c c) , \end{gather}

(4.19)\begin{gather} \frac{\partial a}{\partial t} = D_a(A^0)\left(\frac{\partial^2 a}{\partial z^2} - k^2 a\right) + \frac{\textrm{d} D_{a}(A^0)}{\textrm{d}A^0} \left(2\frac{\partial A^0}{\partial z}\frac{\partial a}{\partial z} + a\frac{\partial^2 A^0}{\partial z^2}\right) \nonumber\\ \hspace{-8pc}\qquad - Da (A_0 b + a B_0) -\psi\frac{\partial A^0}{\partial z} , \end{gather}

(4.20)\begin{gather} \frac{\partial b}{\partial t} = D_b(B^0)\left(\frac{\partial^2 b}{\partial z^2} - k^2 b\right) + \frac{\textrm{d} D_{b}(B^0)}{\textrm{d}B^0} \left(2\frac{\partial B^0}{\partial z} \frac{\partial b}{\partial z} + b\frac{\partial^2 B^0}{\partial z^2}\right)\nonumber\\ \hspace{-8pc}\qquad - Da (A_0 b + a B_0) -\psi\frac{\partial B^0}{\partial z} , \end{gather}

(4.21)\begin{gather} \frac{\partial c}{\partial t} = D_c(C^0)\left(\frac{\partial^2 c}{\partial z^2} - k^2 c\right) + \frac{\textrm{d} D_{c}(C^0)}{\textrm{d}C^0} \left(2\frac{\partial C^0}{\partial z} \frac{\partial c}{\partial z} + c\frac{\partial^2 C^0}{\partial z^2}\right)\nonumber\\ \hspace{-8pc}\qquad + Da (A_0 b + a B_0) -\psi\frac{\partial C^0}{\partial z} , \end{gather}

where $\psi \equiv -Ik{\hat \psi }$, $\varphi \equiv -Ik{\hat \varphi }$. Equations (4.17)–(4.21) should be supplemented by the boundary conditions for disturbances

(4.22a–e)\begin{equation} z ={\pm} H:\quad\psi = 0,\quad \varphi = 0 ,\quad \psi = 0 \quad\frac{\partial a}{\partial z}=0 , \quad\frac{\partial b}{\partial z}=0 ,\quad\frac{\partial c}{\partial z}=0 . \end{equation}

To find stability conditions, the equations for disturbances (4.17)–(4.21) and (4.22a–e) are numerically integrated together with the equations for the base state (3.1)–(3.3) to compute the growth rate $\lambda _L (z,t)$ for a given wavenumber $k$. Repeating the calculations for varying $k$ enables us to describe the growth of small disturbances over a range of wavenumbers. As is known, the growth rate for time-dependent stability problems can be defined in a different way. Based on the fact that the acid in the system is the most active agent, we define the growth rate similar to a Lyapunov exponent

(4.23)\begin{equation} \lambda_L (z_0, t) = \frac{1}{N}\sum^N_{j=1}\frac{1}{\Delta t} \ln\frac{a_{j} (t+\Delta t, z_0)}{a_{j}(t, z_0)}, \end{equation}

where $\Delta t$ is the integration time step and $N$ is the number of independent realizations (typically 10–15). Because $\lambda _L$ is sensitive to the given initial data, each independent integration starts from white noise with an amplitude of less than $10^{-4}$. We fix the occurrence of instability to the time when $\lambda_L (t)$ averaged over $N$ realizations changes sign from negative to positive. The integration time step $\Delta t$ is not constant and changes in accordance with the Courant rule so that the explicit scheme would be stable.

We can apply a linear analysis only to instabilities that develop locally. For example, the SW convection occurs after a global bifurcation, which looks like a catastrophe for the system. Obviously, the study of stability boundaries for infinitesimal perturbations, in this case, does not make sense. Therefore, we can only investigate the development of the CDD and DLC instabilities. Thus, it is convenient to calculate the growth rate of the CDD and DLC disturbances respectively as $\lambda _L (z_{cdd},t)$ and $\lambda _L (z_{dlc},t)$, i.e. at points $z_{cdd}$ and $z_{dlc}$, which coincide with two minima of the density profile (see figure 3).

Figure 7(a) presents the instantaneous growth rates $\lambda _L(z_{cdd},t)$ as a function of the wavenumber $k$ calculated for the CDD convection. This instability develops below the initial contact surface. The minimum of the neutral curve corresponds to the wavenumber $k_{cdd}\approx 4.4$ at time $t_1\approx 0.156$. The figure shows the development of a pure type of instability, which is not affected by other disturbances. This is achieved by considering only those perturbations that are localized inside the potential well. All other disturbances during the calculation have been artificially nullified. Figure 7(b) presents the growth rates $\lambda _L(z_{dlc},t)$ calculated for the pure DLC instability developing above the initial contact surface. This instability starts at $k_{dlc}\approx 0.75$ and $t_2\approx 0.253$.

Figure 7. Real parts of the growth rate $\lambda _L$ of the pure CDD instability (a), the pure DLC instability (b) and the mixed-mode instability (c) are illustrated at different times. The initial concentrations are $\gamma _a=0.667$, $\gamma _b=0.66$.

Although the potential barrier in the form of dense immobile fluid separates the sublayers where two instabilities develop, the instabilities still can influence each other through a diffusion mechanism. To demonstrate the mode interaction, we have recalculated the growth rate curve for both modes, allowing any disturbances to develop in the system. Figure 7(c) presents the growth rates $\lambda _L(z_{cdd},t)$ and $\lambda _L(z_{dlc},t)$ for $t=1$. In this case, the instability region for the CDD convection has significantly expanded due to the longer modes. The DLC mode develops exactly in the same wavelength range. Thus, the influence of the DLC on the CDD instability looks obvious. It is provided by a relatively mobile acid, which easily penetrates the potential barrier and invades the potential well. Thus, it is the acid diffusion that is responsible for the nonlinear mixed-mode instability. To be precise, the DLC instability develops almost independently, since the low-mobility salt and base cannot penetrate the potential barrier separating the sublayers.

5.1. Numerical solution technique

The problem (2.24a,b)–(3.7) is solved numerically by a finite difference method. We consider the Hele-Shaw cell with aspect ratio almost as in the experimental cuvette: $L=20$, $H=37$ (20 per 74 units of dimensionless length). The equations and boundary conditions are approximated on a rectangular uniform grid with $101\times 371$ nodes using a second-order approximation for the spatial coordinates. The one-sided and central differences are used to approximate the time and coordinate derivatives, respectively. The choice of this resolution is based on our previous experience with the numerical simulation of buoyancy-driven flows and is defined by the characteristic size of the arising structures (Bratsun et al. Reference Bratsun, Kostarev, Mizev and Mosheva2015, Reference Bratsun, Mizev, Mosheva and Kostarev2017). The nonlinear equations are solved using an explicit scheme. The Poisson equation is solved by the iterative Liebman successive over-relaxation method at each time step: the accuracy of the solution is fixed to $10^{-4}$. The initial condition for the velocity is determined by a random distribution of the streamfunction field with an amplitude of no more than $10^{-3}$.

Generally, a large Schmidt number allows us to omit the left part of the equation in (2.30). However, in terms of a numerical procedure, it is easier and more efficient to solve a non-stationary motion equation than to iterate while solving the Poisson equation. That is why a nonlinear term in (2.30) is not omitted, although this imposes a special condition on the time step

(5.1)\begin{equation} \Delta t = \frac{\Delta x^2}{2(2+\max(|\varPsi|, |\varPhi|))}. \end{equation}

The following fixed parameters have been used for all the simulations performed in this paper: $Sc = 317$, $Da = 10^3$, $R_{a}=3.18\times 10^5$, $R_{b}=3.82\times 10^5$, $R_{c}=5.1\times 10^5$.

To characterize the nonlinear dynamics resulting from the coupling between chemical reactions and hydrodynamic flows, various types of measurements have been made during the numerical simulations. At successive times, the concentration fields can be spatially averaged along the transverse coordinate $x$ to yield the averaged density profile

(5.2)\begin{equation} \rho_z (z,t) = \frac{1}{L} \int_{0}^{L}\left(A(x,z,t) + \frac{R_b}{R_a} B(x,z,t) + \frac{R_c}{R_a} C(x,z,t)\right) {\textrm{d} x} . \end{equation}

The transverse-averaged profile (5.2) gives information on the speed and intensity of the invasion of the reaction front and convective currents into the lower layer.

If the chemical reactions affect the fluid flows, then convection can also significantly affect the reaction rate, thus ensuring more intensive mixing of the reacting substances. In that respect, an important measurement is provided by the reaction rate computed in terms of the area of the reacted zone (i.e. the number of points where $C(x,z,t)$ is larger than an arbitrary threshold $C_0$) normalized by the area of the system $2HL$ versus time

(5.3)\begin{equation} R(t) = \frac{1}{2HL}\int^{H}_{{-}H}\int^{L}_{0}\theta(C-C_0) \,\textrm{d}z \,{\textrm{d} x} , \end{equation}

where $\theta$ stands for the Heaviside step function. This function returns zero everywhere, except for the areas with a salt concentration $C$ exceeding $C_0$. In the latter case, $\theta = 1$. Typically, the threshold was defined to be $C_0=10^{-3}$.

In what follows, we present the results of direct numerical simulation of buoyancy-driven flows governed by (2.24a,b)–(3.7). We consider here two fundamentally different situations in which reaction–diffusion–convection processes occur under the control of convection ($K_{\rho }<1$) or under the control of diffusion ($K_{\rho }>1$).

5.2. Shock-wave convection

To be specific, we consider the case $\gamma _a=0.667$, $\gamma _b=0.56$ (between points $c$ and $d$ on the stability map shown in figure 4). Mass transfer processes in this area proceed under the control of convection, because the reaction zone is much lighter than the upper layer, and the estimate for the parameter (3.7) gives a value $K_{\rho }=0.97$. Figure 8 shows changes in the integral reaction rate $R(t)$ for this case. One can see that almost complete mixing of salt occurs already after the passage of half a unit of time after the solutions were brought into contact (approximately 150 s). This indicates the extremely intense mass transfer that occurs in the system and the large-scale movement of the fluid throughout the cell.

Figure 8. Variation of the integral reaction rate $R$ as a function of time.

Figure 9 illustrates the evolution of the density profile obtained within a base reaction–diffusion state ($t=0.3$) and in the context of the nonlinear theory ($t=0.7,2,3$). If the linear theory makes it possible to demonstrate the structure of the reaction zone with two minima (see also figure 3), then the nonlinear evolution clearly shows the formation of a SW structure already at the early stage of the reaction–diffusion–convection processes. The pattern is formed as a result of the collapse of two depleted zones of low density. Finally, we note that the SW consists of two density levels with a thin transition region and moves towards a denser lower layer. This figure clearly demonstrates that the analysis of the base-state density profile, as well as the linear stability analysis associated with it, become meaningless. The SW mode is triggered as a result of a global bifurcation. It is interesting to note that the density profile at the bifurcation value $K_{\rho }=1$ can be represented as a homoclinic trajectory (see figure 5c).

Figure 9. Time evolution of a transverse-averaged total density $\rho _z (z,t)$ during the development of the SW convection. The initial concentrations are $\gamma _a=0.667$, $\gamma _b=0.56$. The reaction-induced buoyancy number is $K_{\rho }=0.97$.

Figure 8 clearly illustrates the fact that the system evolution passes through two qualitatively different stages. The first stage is characterized by the exponential growth of the mixing region. In the second stage, the mixing region grows monotonically. Figure 10 shows the density $\rho$ for sequential time moments. The first SW formation stage is presented in (a). At the very beginning of the system evolution, when the convection is still not strongly involved in the mass transfer processes, the densities of the upper and lower layers away from the reaction zone are close to each other and the isopycnal line ($\gamma _a=0.667, \gamma _b=0.56$). Along with that, the reaction–diffusion processes create two layers of reduced densities which are formed from the burning out of the acid and base in the reaction zone, resulting in the accumulation of the heavy salt in the centre of the reaction zone. Since, with the value of the reaction-induced buoyancy number $K_{\rho }<1$, the density of the reaction zone is lower than the density of the upper layer, the entire reaction zone almost immediately after bringing the fluids into contact begins to float up (see figure 10, $t=0.1$). One can see the formation of powerful plumes that emerge under the influence of the buoyancy force (figure 10, $t=0.3$).

Figure 10. Evolution of the density $\rho$ with time showing the development of the SW structure. The frames from left to right and from top to bottom correspond to times $t = 0.1$, $0.2$, $0.3$, $0.4$, $0.5$, $1.0$, $1.5$, $2.5$, $3.0$ and $3.5$, respectively. The integration domain is $0\leqslant x\leqslant 20$, $-37\leqslant z\leqslant 37$. The initial concentrations are $\gamma _a=0.667$, $\gamma _b=0.56$. The reaction-induced buoyancy number is $K_{\rho }=0.97$.

The intense mode described above cannot last for a long time in a closed reactor. The continuous SW development requires a constant input of fresh reactants. Therefore, the evolution of the system becomes calm, which is illustrated in figure 10(b). After approximately $t=0.5$, the upper part of the system becomes quasi-homogeneous, while the lower part, which does not participate in the convective movement, remains homogeneous (figure 10, $t=1$). The system is split into two zones with typical density values divided by the narrow transition zone. A reduced-weight reaction zone continues to be produced, floating up either at the edges of the cell or in certain places. Thus, the fresh acid continues to flow to the contact surface with the lower base layer. Finally, one can observe a SW structure with an almost planar front and nearly discontinuous change in density across the wavefront. This wave propagates fast compared to the characteristic diffusion times and separates the motionless fluid and the area with anomalously intense convective mixing (figure 10, $t=1.5, 2.5, 3, 3.5$).

5.3. CDD convection

If the control parameters are set so that the density of the reaction zone is higher than the density of the upper layer, then a completely different scenario for the development of instability is observed. Figure 11 shows the dynamics of the density field for $K_{\rho }=1.08$ (between points $b$ and $c$ on the stability map shown in figure 4). Note that the system evolution begins with the base reaction–diffusion state, in which the fluid remains at mechanical equilibrium (figure 11, $t=1$). In this state, there is a local maximum located just below the initial contact line $z=0$. Over time, two types of convective perturbation begin to grow in the system.

Figure 11. Evolution of the density $\rho$ with time, showing the development of the CDD convection coupled with the DLC instability. The frames from left to right correspond to times $t = 1$, $2$, $3$, $4$ and $5$, respectively. The integration domain is defined by $0\leqslant x\leqslant 20$, $-37\leqslant z\leqslant 37$. The initial concentrations are $\gamma _a=0.667$, $\gamma _b=0.66$. The reaction-induced buoyancy number is $K_{\rho }=1.08$. The wavelength of the CDD cell and the width of the potential well are determined as shown in the inset, they are used to obtain the theoretical curves in figure 13 of Part 1 (Mizev et al. Reference Mizev, Mosheva and Bratsun2021).

In the upper layer, favourable conditions are created for the initiation of the DLC instability, which usually occurs when two solutions are mixed with a rapidly diffusing solute at the top and slowly diffusing solute at the bottom. In this case, these are acid and salt, respectively.

Below the maximum, a local convective instability arises, which is enclosed inside a potential well. As shown above, the density pocket and conditions for instability development are created under the CDD effect. If this effect is not taken into account, then the liquid in this area remains at rest. That is why this type of instability is called CDD convection (Bratsun et al. Reference Bratsun, Kostarev, Mizev and Mosheva2015). The growth of disturbances leads to the development of a spatially periodic convective structure (figure 11, $t>2$). The integration domain, which corresponds to the width of the experimental cavity, is rather narrow, and therefore only a few cells fit it. However, the numerical simulation with a wider integration domain shows that there exists a periodic system of chemoconvective cells. The obtained pattern is fundamentally different from the usual disordered process of fingering observed commonly in such systems. Most of all, the structure in figure 11 looks like Rayleigh–Bénard convection in a horizontal layer.

An important feature of the system is the independence of the instabilities from each other: they appear in different parts of the medium, are separated by a stably stratified diffusive layer and practically do not interact with each other at the beginning of evolution. Over time, the DLC convection begins to affect the CDD convection by transmitting diffusion signals through the potential barrier separating the two instability zones. Since diffusively coupled instabilities at the beginning have different wavelengths, this causes a LW modulation of one pattern by another. This situation has been studied in detail in Bratsun (Reference Bratsun2019).

The potential well defining the CDD convection zone slowly expands because of the diffusion of the reactants (compare times $t=2$ and $t=5$ in figure 11), and the downward motion of heavy salt-enriched cellular structures follows this expansion. This makes it possible to classify this mass transfer regime as a regime under the control of diffusion.

5.4. Convection cell inside two-layer system

The CDD convection gives a great example of how a localized convective structure can be formed in a medium that generates statically stable potential wells of a density field. It is worth noting that self-sustaining convection localized in the bulk of an almost stationary fluid looks very unusual for researchers who have been educated within the paradigm of thermal convection. In the latter case, the medium, as a rule, is heated from the boundary side and convective motion quickly expands to the entire accessible region of the cuvette. The chemical reaction that destabilizes a liquid medium is fundamentally a local tool. By using the CDD instability, we can create, at least in numerical simulation, a convective cell that is statically stable and works autonomously for a sufficiently long time.

To illustrate the above statement, consider a two-layer miscible system placed in a square Hele-Shaw cell (40 per 40 units of dimensionless length), where the lower layer is an aqueous solution of the base, and the upper layer is an aqueous solution of a neutral solute that does not react with alkali. The density of the upper layer is less than the density of the upper. Thus, at the very beginning, the system is stable to perturbations of the RT instability. Suppose that the acid diffuses into the system from above in a narrow band defined by $x=[23,33]$. Then, the reaction occurs only at the intersection of this band and the contact interface. Further, the diffusion of reactants results in the bulk reaction. Figure 12 shows the vorticity contours and total density field at the time $t=3$. One can see that edge effects affect the shape of the reaction front. The CDD instability appears below the potential barrier, while two rising plumes of the DLC convection are developed above the barrier. The vorticity field shows that the fluid motion is localized in a square convective cell measuring 10 by 10 units. The fluid around almost does not take part in the movement.

Figure 12. Vorticity field $\varPhi$ (a) and the total density $\rho$ (b) show the occurrence of a localized convection cell inside a two-layer system when the acid solution enters the system from above in a narrow band defined by $x=[23,33]$. The integration domain is defined by $0\leqslant x\leqslant 40$, $-20\leqslant z\leqslant 20$. The initial concentrations are $\gamma _a=0.667$, $\gamma _b=0.66$. The reaction-induced buoyancy number is $K_{\rho }=1.08$.

5.5. Comparison between experimental and theoretical studies

By inspecting the experimental data presented in Part 1 and the theoretical results presented in Part 2, we can conclude that there is a good agreement between experimental observations and theoretical modelling. The principal stability maps obtained experimentally (see figures 6(b) and 14(a) in Part 1 (Mizev et al. Reference Mizev, Mosheva and Bratsun2021)) and theoretically (figure 3) coincide qualitatively and quantitatively. The modes of chemoconvection observed experimentally (see figures 4 and 12 in Part 1 (Mizev et al. Reference Mizev, Mosheva and Bratsun2021)), are almost indistinguishable from the regimes obtained in numerical simulation (figures 10 and 11).

Figure 13 demonstrates important quantitative comparisons between experimental data and numerical simulations. Figure 13(a) shows the change in the width of the potential well $d$ over time. As noted in Part 1 (Mizev et al. Reference Mizev, Mosheva and Bratsun2021), the method of Fizeau interferometry makes it possible to qualitatively visualize the spatial distribution of species concentration. This allows comparison with theory. One can see from the figure that, except for the very initial stage, the width of the well increases linearly with time. Good agreement between theory and experiment indicates that the theoretical reaction–diffusion model works correctly.

Figure 13. Comparison of experimental data and results of numerical analysis: the growth of the potential well width $d$ over time (a); the variation of the wavelength of chemoconvective cells $\lambda$ scaled by the potential well width with time (b); the speed of the reaction front $V$ in the mode controlled by diffusion ($K_{\rho }>1$) and convection ($K_{\rho }<1$) (c).

Figure 13(b) shows the evolution of the wavelength $\lambda$ of chemoconvective cells scaled by the potential well width $d$. One can see that the cells appear earlier in the experiment than in numerical simulation. This may be due to the fixed amplitude of the perturbations of the streamfunction that we used as the initial state of the system (less than $10^{-4}$) and the slow growth of perturbations at the very beginning of evolution. It is interesting to note that the wavelength of convective cells grows according to the root law $t^{-1/2}$, while the width of the potential well, where the cells develop, grows linearly (figure 13a).

The important characteristic of a shock wave is its velocity. The system enters the SW mode only if the wave acquires a speed higher than some critical value $c^*$. We have previously shown in Bratsun et al. (Reference Bratsun, Mizev, Mosheva and Kostarev2017) that the critical velocity can be defined as

(5.4)\begin{equation} c^{\ast} = \sqrt{Sc}=\sqrt{\frac{\nu}{D_{a0}}} . \end{equation}

In our problem, $c^*$ plays the same role as the speed of sound in gas dynamics. This means that, if the density wave moves faster than $\sqrt {Sc}$, it behaves as a subsonic analogue of the shock wave in gas. The velocity (5.4) estimated for our problem as $c^*\approx 17.8$ ($0.05\ \textrm {mm}\ \textrm {s}^{-1}$ in dimensional units) is in perfect agreement with experimental data. We found that, as soon as the propagation velocity of the SW falls below this value, the wave immediately stops and is replaced by a common fingering under the DC mode. Figure 13(c) shows a comparison between the reaction front velocity at the very initial moment of evolution depending on the value of $K_{\rho }$. One can notice that the experimental values slightly exceed the velocity values obtained in the numerical analysis. This difference may be related to the value of the Damköhler number, which was fixed in the simulations. Since we do not know the exact value of the reaction rate $K$ in the medium under consideration, we fixed it to the value $10^3$. Most likely, in a real system, this value is a bit higher and the reaction in the experiment proceeds more vigorously.