Convective instabilities derived from dissipation of chemical energy

Pure reaction-diffusion structures and propagation have been extensively studied and are well-known from the basis of a Turing instability.^1,2 The presence of a reaction front is also a well-known feature in many physical, biological, and chemical processes. For example, a simple chemical reaction type process of A + B → C will exhibit a front, a spatially-localized region where production of C can be achieved, provided diffusing reagents are initially separated in space.³ Even if A and B are initially not separated, a lateral instability can still emerge if the prevailing reaction is excitable, e.g., autocatalytic.⁴ When a spatially extended system becomes unstable, the resulting dynamics have sensitive dependence on whether the system is convectively unstable, in which case perturbations grow in time but convected away fast enough that they die at each fixed position in the lab frame considered, or absolutely unstable, in which there arises a perturbation and a location where this perturbation does not decay. With an autocatalytic reaction system, a perturbation using an autocatalyst does not decay and will grow until diffusion or convection-limited. The first stage in understanding pattern formation in such reactive flows is to understand the description and calculation of the properties of the reaction zone. This involves answering the question of where and at what rate does the reaction product appear. This problem has been studied with respect to Liesegang-ring band formation in which simplified theories have predicted band formation solely from properties of the front and phenomenological theories of diffusion-limited aggregation.^5,6 The width of the reaction zone determines gradients of several physical and physicochemical parameters and hence the relevant driving forces of wave and pattern generation. Simple bimolecular kinetics have been used for reaction dynamics, which is too simplified a model and is incorrect for autocatalytic exothermic fronts.⁷ Numerical simulations have, thus, not been accurate in predicting the structure, especially in the wake of the propagating front.

Convective flows due to thermocapillary Marangoni-type convection coupled to Rayleigh-Taylor, Rayleigh-Benard, or double-diffusive types of hydrodynamics flows can be triggered by spatial gradients of concentrations and temperature.^8,9 Such behavior has been observed in real-life phenomena in engineering systems, geological formations, supernovae dynamics, CO₂ sequestration, earth mantle convection, and atmospheric and ocean layering convective motions.¹⁰ 

This article assesses the active role that chemical reactions can have on parameters of flow dynamics in chemical reactive flows. We have, previously, been able to detail the mechanism of the “drive” reaction that drives the hydrodynamic instabilities reported here.^11–13 The reaction used is the autocatalytic chlorite-thiourea reaction.^11,14 Autocatalytic reactions coupled to diffusion can generate traveling fronts of constant speed whereby products of the reaction invade and progressively consume reactants in space and time.¹⁵ In the absence of a gelled environment, convection can deform such fronts due to density differences between reactants and products. Density differences can result in differences in molar volumes of products and reactants. This is especially so in the chlorite-thiourea reaction where products such as sulfate and chloride are dense ionic products while reactants are lighter organic reagents. The exothermicity of the reaction, however, can temporarily “flip” this density relationship until thermal equilibrium is attained.¹⁶ Thus, convective motions will appear around descending fronts since we have a dense solution overlying a lighter one (reactants). This is an unstable front since any deformity of the surface will tend to grow. A Rayleigh-Taylor instability (RTI) in which denser solution sinks into the rising less dense reactant solution will deform the interface into density fingering since the density gradient points opposite to the body force. This is due to the buoyancy-induced convective deformation of a reaction-diffusion front. The combination of convection to a reaction-diffusion front results in reaction-diffusion-convection dynamics. Our preliminary experiments have shown that addition of heat and, subsequently, convection (both Marangoni and Rayleigh-type) has not produced more complexity, but more cooperativity.^17–19 This is counterintuitive to what we would have expected.²⁰ This fascinating cooperativity has resulted in hitherto previously unobserved spatiotemporal inhomogeneities, Rayleigh-Taylor²¹ and Richtmeyer-Meshkov instabilities,²² highly ordered and periodic thermal plumes,¹⁹ Turing-like concentric patterning, and an independent series of convective tori with well-defined wavelengths, all derived from a single batch reaction.¹⁸ 

Chlorite-based oxidation reactions show the greatest complexities due to, among others, interfacial effects which can introduce stochasticities.²³ Universally, after an induction period, all thiocarbamide-oxyhalogen systems produce a large pH change (as H⁺ ions are rapidly produced) as well as large heat generation. The initial oxidation of the S(-II) species to S(0) (e.g., sulfenic acids) has a negative acid feedback loop, but subsequent oxidations all the way to S(VI) produce a positive acid feedback loop.¹² The only prerequisite for the observations of the full range of exotic dynamics is that the oxidant (chlorite) should be in stoichiometric excess over the reductant.¹¹ The advantage, then, is that since [H⁺] is allowed to vary and is determined by the reaction's stoichiometry, then an acid propagation wave can be observed and used to follow the lateral instability through an acid-base indicator. This approach has also been utilized in following waves generated in the chlorite-tetrathionate (CT) reaction.^24–26 The stoichiometry of the $ClO2−$-thiourea reaction is^11,12,27

We have derived a general mechanism to account for most chlorite-based oxidations of sulfur compounds, and it appears the chlorite-thiourea reaction can be explained, in general terms by this model. Reaction (R2) is a general stoichiometric equation for the oxidation of thiocarbamide compounds by chlorite.¹¹ The oxidation products are sulfate and a urea-type residue, R₁R₂$C=O$⁠,

These reactions show common clock reaction characteristics in which initially there is a quiescent period with no activity in the reaction indicators. At the end of the induction period, there is a sudden production of sulfate, $SO42−$⁠, chlorine dioxide, ClO₂, and acid, H⁺. The acid and sulfate are formed simultaneously from the oxidation of the sulfonic acid intermediate, $R1R2′CSO3H,$

Chlorine dioxide is formed by an extraneous reaction between the excess chlorite, $ClO2−$⁠, and left-over autocatalyst: hypochlorous acid, HOCl,²³ 

HOCl, being the autocatalytic species, controls the rates of reaction and wave propagation under isothermal conditions. The mechanism for autocatalysis involves the reactive intermediate Cl₂O₂ leading to quadratic autocatalysis²⁸ 

followed by

So, for almost all our modeling studies, the rate of reaction can be simplified into Rate = k[A][B] for a standard A + B reaction where B represents the autocatalyst.

Apart from H⁺, $SO42−$⁠, and ClO₂, the most important product parameter in terms of wave propagation dynamics is a large and almost instantaneous heat evolution. The heat generated is very high, such that the reaction of thiourea and chlorite has a reaction enthalpy in excess of –1170 kJ mol⁻¹.^14,29 The transformation from purely organic reagents to mostly inorganic products; sulfate, H⁺, Cl^– results in high energies of solvation and an increased solution density. By adjusting initial reagent concentrations, temperature jump for our reaction system could be maintained between ΔT = 2.2–3.8 K. This is the range that produced the most exotic wave and front dynamics. The unreacted solution far ahead of the wave front maintained the room temperature while the leading edge of the wave front could be 2.2–3.8° above this temperature. The reaction zone is a thin sliver of approximately 2 mm.³⁰ The experimentally-determined reaction zone and temperature profile are shown in Figs. 1 and 2. This was obtained using thermocouples lined along the wave front.

There are several physical parameters that can be utilized to follow the development of patterns and/or waves in this reaction system. Apart from neutral density polystyrene beads, there are four possible indicators that can be used to follow the traveling wave and the subsequent spatiotemporal patterns: acid-base, BaCl₂, ClO₂ color, and freshly-prepared starch. Each has a specific role that the others cannot duplicate. BaCl₂ gives a white precipitate of BaSO₄ in the presence of sulfate, which is produced in reaction (R3).³¹ This article dwells only on the use of acid-base indicators. This method is profiled below: in unbuffered solutions, protons produced by the reaction can lead to large changes in pH of the solution. Typical pH values are in the range of approximately 5.5 in reactant solutions and 1.5–3.5 in product solutions [initial reactant concentrations determine the amount of H⁺ ions formed, see stoichiometry (R2)]. Acid-base indicators are best suited for the study of front propagation and spatiotemporal patterns due to their vivid color at low micromolar concentrations. Acid-base indicators are also very good at being able to characterize the reaction zone of the wave front: intensity and degree of color change with indicators methyl red and bromophenol blue is indicative of the extent of reaction. This gives a very good estimate of the reaction front (which we have so far estimated at about 2–3 mm;¹⁸ exact values for conditions utilized for our experiments were obtained in this study, see Fig. 1). Acid-base indicators are not effective in aiding in the observation of hydrodynamic patterns since they form a homogeneous solution with a reaction solution.

The following analytical grade reagents were used without further purification: thiourea (Aldrich); perchloric acid 70%–72%, bromophenol blue, methyl red, and methyl orange (Fisher). The structures of these three dyes used are shown below. Sodium chlorite is sold in its technical grade form (about 80%). A single recrystallization (ethanol-water mixture) brought the assay to >98%. The analysis of sodium chlorite was performed iodometrically by adding excess acidified iodide and titrating the liberated iodine against standard thiosulfate with freshly-prepared starch as indicator.³² 

Stock solutions of sodium chlorite were prepared fresh for each set of experiments and stabilized with 0.001M sodium hydroxide if the experiments were going to be run in high acid environments. The sodium chlorite solutions were also stored in dark Winchester vessels wrapped in aluminum foil to reduce decomposition by light. The observation vessel used was a modified Hele-Shaw cell 260 mm long and a thickness of either 5 or 10 mm machined out of plexiglass. Quartz vessels gave stochastic triggering and thus quartz vessels were not utilized. The total depth of the vessel was 30 mm, but reaction solutions were filled to not more than 20 mm in depth, leaving a 10 mm head. This decreased the degree of possible evaporative cooling that could affect thermocapillary convection.

The chlorite, thiourea, acid, and/or indicator were thoroughly mixed before being poured into the reaction vessel. Chlorite was added last by a rapid delivery pipet. The wave could self-initiate, but the induction times for clocking were stochastic. The time taken before clocking was very important in determining wave dynamics, especially velocity and shape. For example, solutions that clocked early gave a faster-moving wave front from the point in initial perturbation. After solutions were thoroughly mixed, they were allowed to settle until all physical ripples disappeared and were triggered by the addition of a 0.025 ml micropipette drop of a solution containing HOCl. Due to autocatalytic acid generation at the wave front, pH indicators are ideally suited to study the wave front dynamics and specifically the reaction zone of the wave front. The vivid color change of these indicators can be easily enhanced and quantified using MATLAB. This also allows for very precise measurements to be made regarding the width and velocity of the wave front. By utilizing a unique combination of indicators, methyl orange (pKa = 3.4) and bromophenol blue (pKa = 3.8), we are able to visualize the exotic structures that are generated from the interaction of thermocapillary and thermogravitational forces as well as the standard reaction-diffusion propagation mechanisms.

Previously, we had established that, in shallow layers, 1.7 mm<d < 2.3 mm, wave front propagation is proportional to the temperature jump at the wave front, ΔT. The temperature jump is proportional to the heat production by the reaction Q. Q, in turn, is a function of the extent of the reaction, ξ,

ΔH_R is the enthalpy change of the reaction and c_i0 represents initial concentrations subject to stoichiometric constraints that the ratio of oxidant to reductant is greater than 2. Thus, at constant initial concentrations, the variable which controls amount of heat generated within the reaction zone is the extent of reaction, ξ. ξ varies from zero just in front of the wave to 1.0 at the end of the reaction zone. The maximum value of ΔT can be derived from the heat capacity of the reaction solution in the limit of extent of reaction approaching unity. The highest temperature is recorded just behind the wave front, and, subsequently, a drop in temperature is recorded based on Newton's law of cooling. In shallow layers, the wave accelerates until it attains a constant lateral velocity. Since the Marangoni effect dominates in shallow layers, the wave velocity is, thus, determined by the rate of reaction and thus reaction kinetics. This bistable autocatalytic chlorite-thiourea reaction has the extent of reaction that is a nonlinear function of the instantaneous chemical reactant concentrations. The autocatalysis term will dominate the kinetics of the reaction; this is the R5 + R6 sequence. The reaction that delivers acid that determines the position of the front is the oxidation of the sulfonic acid in reaction (R3),

Thus, the rate of formation of acid is dependent on the rate of formation of HOCl. The R5 + R6 can be summarized as

The rate of formation of HOCl is

Equation (2) can be integrated as an initial value problem, giving the concentration of HOCl at any time as

Equation (3) is a sigmoidal trace which is limited by the diffusion of reactants in the reaction zone and the depletion of the substrate. Thus, we expect an ever-increasing rate of reaction which tails off when thiourea is depleted. These kinetics are much more complex than those expected from the standard A + B → products, which has been the staple of most previous studies of chemical reactive flows.³ 

Our studies have three different waves in one wave front. The leading front is the acid wave, determined by the change in pH indicator. This is derived from reaction (R3). The second wave can be visually observed by barium sulfate precipitation in reaction solutions spiked with barium chloride. This is also derived from reaction (R3), but the proton's diffusion coefficient means that the BaSO4 wave trails behind the reaction front. The third wave is of yellow chlorine dioxide, derived from reaction (R4). The chlorine dioxide wave could be enhanced by the addition of freshly-prepared starch indicator preserved with HgI₂. The acid front, determined by the acid-base indicator, was then annihilated by bleaching from chlorine dioxide. Thus, a colorless solution emerges way behind the leading wave front.

There are several parameters that are important in determining observed spatial inhomogeneities, most important being the initial reagent concentrations as they determine the amount of heat produced which will fuel wave propagation. Our initial studies of wave propagation had involved the following initial concentrations: [NaClO₂]₀= 4.00 × 10⁻³M; [CS(NH₂]₀ = 2.00 × 10⁻⁴M. Upon triggering the reaction, the temperature increases, with ΔT = 3.3 K. This was also collaborated through calculations from the expected heat of reaction, heat capacity of the vessel, and assuming reaction solution has the same specific heat as water. The solution depth was varied from 1.7 mm to 3.5 mm. Depths shallower than 1.1 mm could not deliver a completely flat surface due to the surface tension between the plexiglass and the reaction solution. Depths of 1.3 mm <d < 2.1 mm delivered a lateral instability of an acid with no structure. Wave velocity commenced with a rapid acceleration followed by a steady wave velocity as would be predicted by the rate of formation of HOCl in Eq. (3). The isothermal density change was measured as Δρ = 2.8 × 10⁻³ g cm⁻³. The change in surface tension with temperature, as measured by the capillary rise method, ∂σ/∂T = −9.5 × 10⁻⁵ Nm⁻¹ K⁻¹. A negative value for ∂σ/∂T is expected for most liquids. ∂σ/∂T was measured only for the product solution. The expansion coefficient for the product solution, ∂ρ/∂T, was determined experimentally as −0.4 kg m⁻³ K⁻¹.

In previous work performed in our laboratory, the pH value of solution varied from 5.5 to approximately 2.6. Such a range was amenable to the observation of acid front by the use of methyl orange and methyl red indicators. In the reported experiments in this article, the final pH varied between 2.6 and 3.8. While methyl orange/red could catch the lowest parts of the pH range, they were not as effective in the higher ranges. Thus, a mixture of the indicators, methyl orange and bromophenol blue, was used to determine not only the wave front but also the extent of the reaction in the reaction zone before bleaching by chlorine dioxide occurred. Based on the mixture of pH indicators used, the unreacted solution had a bluish-green color, the wave front is red, and it varies in reddish intensity based on the extent of reaction. Chlorine dioxide formation in the wake of the acid front, as mentioned previously, bleaches the indicators to a white color.

All experiments reported here were run with the following initial reagent concentrations: [CS(NH₂)₂]₀ = 1.68 × 10⁻⁴ M, [$ClO2−$]₀ = 4.00 × 10⁻³M for an oxidant to reductant ratio of 24:1. This ensured the total consumption of the thiourea. Product solutions contained chlorite (excess), chloride, chlorine dioxide, sulfate, and acid. These experiments were not buffered so that acid variations could be used to determine reaction front and rate. Initiation was solely through the addition of autocatalyst at the edge of the Hele-Shaw cell. Same wave dynamics were observed from electrochemical initiation when a platinum electrode is placed on one end of the vessel and a voltage of 3.0 V is passed through. The autocatalyst, HOCl, was added in the form of a basic solution of chlorine water through a micropipette.

Added at the edge of the Hele-Shaw cell, the autocatalyst immediately commenced the reaction, and this was observed through the immediate acid formation [reaction (R3)] and yellow chlorine dioxide [reaction (R4)]. Due to the positive isothermal density change, the acid wave descended down the cell to the bottom of the cell. The proximity to the edge of the vessel ensured efficient heat exchange of the warm front with the ambient temperature. The efficient heat exchange with the unreacted cold solution ensured its rapid fingering down the cell. This was the initiation of the rest of the wave front dynamics reported in this article. Figure 4 shows a sketch of all the relevant and expected forces in this reacting-diffusing-convecting medium. Our previous studies had shown a dominant Marangoni effect at the surface in shallow layers.

The surface tension gradient on the free surface will induce a top layer which has to be materially balanced by a return back flow at the bottom of the fluid. Villiers and Platten^33,34 modeled this flow for high Prandtl number fluids, Pr ≥ 4, and evaluated a horizontal velocity profile of

where V_x is the horizontal velocity and h is the normalized solution depth. The sketch of this flow is shown in Fig. 3 for depths between 1.7 mm and 3.3 mm. These shallow layers had a dominance of thermocapillary effects. Figure 4 shows fluid flow with a reaction solution with reagent concentrations that gave a front temperature jump of 3.3 K. The thin sliver of product solution pulled forward due to thermocapillary effects is brought to thermal equilibrium with unreacted bulk solution resulting in fingering regime of double-diffusive convection. As the finger proceeds down the bulk unreacted solution, it is swept backward in accordance to expected fluid flow depiction in Fig. 3. Depths and values of ΔT could be adjusted to deliver various patterns and fluid flow dynamics. These patterns include thermal plumes and convective tori with entrainment of fluid flow in concentric circular patterns.

This article reports mainly of fluid flow patterns generated from 20 mm depths of a 260 mm long Hele-Shaw cell. The reaction solution used, coupled with the heat capacity of the apparatus, gave a traveling wave with a maximum temperature jump after the complete reaction, of 2.1 K. This value of ΔT is lower than the temperature jump for solutions that gave convective tori and thermal plumes. Thus, these experiments were performed at higher Eötvös numbers than those previously reported,

This Eötvös number is effectively a ratio of thermogravitational and thermocapillary forces. The lower reduced Marangoni effect delivered hitherto unobserved patterns of self-organization: spatiotemporal inhomogeneities and Rayleigh-Taylor instabilities. Figure 6 shows the first waves generated after the initial initiation of the wave and full descent of the initial finger on the side of the container, x = 0. Finger A (faint) originates from the top of the solution and immediately fingers down the bulk of the solution. As opposed to Fig. 5, the initial wave at the top of the solution immediately starts to sink into the bulk solution without extensive movement in the x-direction. Finger A's forward velocity stalls at 0.858 mm min⁻¹ as the mass of products behind the wave front cools and increase their specific gravity relative to the unreacted bulk solution. The descending wave velocity is oscillatory as it descends through the 4 images [Figs. 6(a)–6(d)]. It commences at 3.78 mm min⁻¹, accelerates to 3.87 mm min⁻¹, and finally decelerates to 3.83 mm min⁻¹ as it approaches the bottom of the vessel. This deceleration is due to the proximity of the other hot region at the bottom of the vessel, and, thus, heat loss and subsequent decrease in specific gravity of the product solution is reduced. The wave at the bottom of the vessel proceeds in the x-direction at a constant velocity of 4.02 mm min⁻¹. The rising plume is derived from the increase in temperature due to the autocatalytic reaction, causing the hot products to rise rapidly through the reagent solution as a volcanic plume. The plume rises through at 3.05 mm min⁻¹ in Fig. 6(b) and slows down to 1.54 mm min⁻¹ in Fig. 6(d) when it spreads horizontally back into finger A. From the expected dynamics of the wave front, we expect the descending finger to be unstable (any defect at the wave front will be magnified) and the ascending plume to be stable (any deformities at the interface will be suppressed).

Figures 7(a)–7(d) show the further development of patterning from Fig. 6. Here, we concentrate on the development of a new finger. For a short period, wave moves with v_x = 2.83 mm min⁻¹. Again, due to the lower ΔT at the wave front, the leading front of the wave quickly dips below the surface of the liquid and is propelled more by reaction-diffusion-convection mechanisms and not as much from the Marangoni instability. The formation of the fingering regime effectively halts the horizontal lateral instability of the effectively hydrothermal wave. Once a minor defect is created at the wave front interface, the denser products overcome the Marangoni convection and cause a rapid descent as the products force the wave front downward. This initial descent commences at 12.1 mm/min in the first one third of the reaction depth. Following this rapid descent, the wave front begins to curl back to the wave front at 1.16 mm/min due to availability of reagents and expected back flow, though not as pronounced as would have been expected from Fig. 4. After reaching the bottom of the container the wave front begins to rapidly propagate in both x-directions moving to the left at 4.70 mm/min and to the right at 2.93 mm/min. Propagation to the left is through a combination of RDC mechanisms and expected back flow (see Fig. 4), while propagation in the positive x-direction is solely through RDC mechanisms.

Figure 8 shows further progression of the hydrodynamical instabilities from Figs. 6 and 7. This is after three iterations of fingering and plume. The patterning changes in chemical wavelength, wave front speeds, and velocities of the plumes and fingers. There are two reactions proceeding in the reaction mixture: the catalyzed and noncatalyzed reactions. Thus, the chemical composition of reactants keeps changing with time. Figure 8 shows the usual initial lateral forward movement near the top of the solution, labeled as “A.” The finger “B” slants in the positive x-direction because of interaction with the plume movement just below it. Thus, the continuation of the plume in “C” has to proceed in the positive x-direction. Near the top of the solution, plume “C” bifurcates into “D” in the negative x-direction and “E” in the positive x-direction. The velocity of lateral instability “D” is much less than that of its positive counterpart “E.” D is slower because of the lower thermal gradient due to incoming instability A. Separation of A and D is maintained for extended periods. Instead of coalescing, they, instead, both finger into the unreacted region created by B and C resulting in very exotic heterogeneities in acid, sulfate, chlorine dioxide, and heat. After lateral instability denoted as “E,” the general sequence shown in Figs. 6 and 7 is repeated of finger and plume alternating. Our observation vessel is 260 mm long and thus not too many of these sequences are observed. Also, as wave approaches the end of the Hele-Shaw cell, boundary and eddy effects from the other wall begin to affect wave movement. In general, wave movement, propagation, and spatiotemporal organization are heavily dependent on vessel size and geometry. Figures 7 and 8 show entrainment of unreacted reagents, surrounded by zones of reacted products, specifically the separation of fingers A and D in Fig. 8 and the isolated reagents in Fig. 6(d). If observed along the z-axis; the system will display spatiotemporal inhomogeneities in the x-axis. If the reaction vessel is wide enough, these spatiotemporal inhomogeneities will cover the x-y plane.

Chemical waves generated in an exothermic or endothermic reaction produce density inhomogeneities in the reaction medium. In considering double-diffusive flow between two parallel plates that are perpendicular to the direction of gravity, the Oberbeck-Boussinesq approximation is used in which fluid density is assumed to depend linearly on two scalar fields, namely, the temperature, T, and the respective concentration of the different solution constituents, c_i. In our nonlinear reaction system, the density will, thus, vary depending upon the extent of reaction and the temperature. The density, ρ, can then be derived from

where

with subscripts “0” denoting initial conditions. α is the thermal expansion coefficient and β is the solutal expansion coefficient. The first term in Eq. (6) represents the temperature dependence of density and the second term the concentration dependence. These two terms, in the system under study, deliver opposing effects, which give rise to the observed multicomponent convection (fingering regime). At isothermal conditions, the first term vanishes, leaving only the second term which gives, in this reaction system, ρ > ρ₀. This inequality will dictate the effect of buoyancy at thermal equilibrium. However, for both reactants and products, the buoyancy factors can be estimated by the Grashof number,

where d is the depth of the solution, κ is the thermal diffusivity, and ν the kinematic viscosity. Since in our experiments we evaluated δρ/δT = −0.4 kg m⁻³ K⁻¹, then when T > T₀, the first term in Eq. (4) will pull the value of density, ρ, to lighter values. With d raised to the third power, buoyancy factors are much more pronounced in deeper fluid solutions.

At isothermal conditions, the product solution has a slightly lower surface tension than the reactant solution. Thus, both reaction-diffusion and thermocapillary convection act in the same direction with respect to wave propagation,

Surface tension, in general, decreases with temperature, and our reaction solutions gave τ = −9.5 × 10⁻⁵ N m⁻¹ K⁻¹. Thermocapillary forces can be estimated from the Marangoni number, Ma[18],

where L is the length of the free surface. Although ΔT changes over a reaction zone of about 2 mm (see Fig. 1), a reasonable Marangoni number can be calculated from our experimental data. The Eötvös number is effectively the ratio of thermogravitational and thermocapillary forces,

where Ra, the Rayleigh number = Gr ⋅ ν/κ is a better measure of the relative strengths of the forces responsible for the formation of the convective torus. In general, we calculate that the Eötvös number has a temperature-independent value of 0.40 for the same vessel and same conditions.³⁴ This means that, for the same temperature gradient, the Marangoni effect should be 2.5 times larger than the thermogravitational-based Rayleigh effects. We expect, theoretically, a strong dominance of thermocapillary effects, with ΔT, over thermogravitational effects, and this was also observed experimentally in a previous publication from our laboratories. It is, thus, easy to investigate our pseudo-two-dimensional thermocapillary-buoyancy flows. Numerical analysis of the strongly-convective thermocapillary flows has shown that a finite temperature gradient remains over the entire free surface, leading to a global flow structure. Theoretical work by Smith and Davis also concluded that thermocapillary convection dominates in such cavity problems in the limit of fixed E₀ and as Ma → infinity.

Our system is a modified form of a Benard-Marangoni (BM) instability with internal heat generation. Convective instabilities driven by either buoyancy or thermocapillary effects have been the subject of numerous theoretical and experimental studies since the initial studies of Lord Rayleigh and Pearson. Considerably less work has been done on the more general and relevant case in which both mechanisms are acting simultaneously. Pioneering theoretical work on this BM instability was performed by Villiers and Platten in which they used acetone with lateral heating.³³ They obtained qualitative agreement with experimental data with respect to the various convective regimes: steady or oscillatory. Velocity profiles were those calculated and shown in Fig. 9. Prandtl number for acetone was fixed at 4.24, and the major bifurcation parameter was the temperature jump imposed at the wave front. They modeled transient convection in a finite rectangular enclosure by using the stream function vorticity formulation of the two-dimensional Navier-Stokes and energy equations and assuming the Boussinesq approximation. By utilizing the relation

and the Marangoni number [Eq. (10)] and the force for a nondeformable flat surface, they could utilize the no-slip adherence boundary condition and the Levich boundary condition

to determine the commencement and nature of the lateral instability. Our reaction-diffusion-convection system requires several other parameters that come into play. Specifically, for a full modeling of the system, one needs evaluation and determination of the following parameters: Atwood number, Lewis number, Damköhler number, Schmidt number, Peclet number, the Nusselt number, and the Soret effect.³⁵ Many phenomena and parameters like the Atwood numbers and Soret effect can be evaluated for our system.^36–38 We determined an Atwood number, A_t, of 1.4 × 10⁻³. In Rayleigh-Taylor instability (RTI), the penetration distance of heavy fluid bubbles into the light fluid is a function of acceleration time scale, Agt,² where g is the gravitational acceleration and t is the time.³⁶ In our system, both natural and buoyant convections are important [see Eq. (4)]. There is no absolute cut-off value, but, generally, with A_t < 0.1, we can assume that the velocity field is solenoidal and thus use the Boussinesq approximation.³⁹ At small Atwood numbers, the development of Rayleigh-Taylor instabilities proceeds initially through an exponential growth of infinitesimal perturbations that correspond to linear stability analysis. At amplitudes of about half of the wavelengths, the RTI saturates, and longer wavelengths take over. There will finally be nonlinear mode interactions for final successive wavelength saturation. We can define changes in temperature and concentration at the front as

T₁ is calculated from the extent of reaction and assumption that there are no heat losses. The Soret effect quantifies the fact that a flux can be generated solely by a temperature difference.^40,41 This effect does not exist in pure fluids, and thus we can define a factor D_T for c₀(1 – c₀). D_T remains temperature dependent exactly as the original D, before the temperature jump. Total flux is the standard Fick's term, combined with a mass separation due to the temperature gradient. These two terms are of opposite signs, and when we are in steady state conditions, the terms will be equal with zero total flux, implying

We define the Soret coefficient as

This effect can be positive or negative, depending on the sign of D_T. We have attempted, initially, to estimate our Soret coefficient. The reaction commences with mainly organic compounds and has inorganic ionic products. We expect an absolute Soret coefficient for organics of approximately 1.00 × 10⁻³ to 10⁻² K⁻¹. In our system, with a ΔT of 1.5–2.5 K, we expect a negligible Soret effect since this will normally approximate to the extent of reaction within the linear region of temperature jump at reaction front. The Soret effect would be sizable, however, in fluid physics of crude oil reservoirs. We also expect thermal conductivity and kinetics constants to be essentially invariant after the temperature jump.

We assumed no convection in the y-axis and accepted a 2D model as plausible for this system with the incompressible velocity field u = (u_x, u_z). The system has length L_x and height/width L_z. Gravity points along the z-axis. We then apply Darcy's law [Eq. (17)] and continuity equation [Eq. (18)] and combine this with the reaction-diffusion-convection evolution equations for the autocatalytic production of the product species, c, and for the elevation of the temperature, T, derived from the exothermic reaction [Eqs. (19) and (20)]

In Darcy's law, density depends on concentrations and temperature as quantified in Eq. (6). “D” in Eq. (15) is the diffusion coefficient. In Eq. (20), α is the Newton heat loss coefficient. Since we worked only in the reaction front, this last term in Eq. (20) was ignored. We need to obtain dimensionless equations and parameters. We proceed to define the following:

The Lewis number is the ratio of thermal diffusivity to mass diffusivity. Le = λ/ρDC_p; and it is also a ratio of the Schmidt to the Prandtl number. Equations (17), (19), and (20) can then be modified to nondimensional forms after removing the primes. Of note are the two forms of Rayleigh numbers, one for solutal buoyancy, R_S, and the other for thermal buoyancy, R_T,

In Eq. (21), sign of i depends on whether wave motion is buoyantly stable or not. Most of the relevant parameters for our system can be calculated and/or evaluated. We assumed the same solution viscosity as that of pure water of 1.00 centiPoise. We evaluated an approximate R_T value of 28.5 (dimensionless). Our limited model managed to model the initial lateral instability velocity. It also predicted fingering instabilities. This is, effectively, the flows shown in Fig. 4. For plumes, boundary conditions were difficult to determine. Several parameters could not be accessed in our system. The experimental procedures involved an initiation at a point source, thus precluding the RTI instabilities that have been observed in other autocatalytic systems such as the chlorite-tetrathionate system. We were interested in the development of a lateral instability as well as the fingering and pluming regimes. Nonlinear chemical kinetics and the limited Hele-Shaw dimensions quickly introduced boundary conditions. By assuming Dirichlet boundary conditions, where the trajectory of the plume is assumed at the bottom of the vessel, the oscillatory nature of the plume speed could be qualitatively but not quantitatively modeled.

The remarkable aspect of our studies is the similarity they may have in stellar evolution.⁴² The understanding of the structure and evolution of stars is a long-standing problem in astrophysics. The central theme is the description of stellar turbulent convection and how it transports heat and energy.

In massive stars, there exists a semiconvection zone, a chemically inhomogeneous zone between the convective core and the radiative envelope.⁴³ The governing equations mirror the equations we have derived here in this article. As the star evolves, a gradient of chemical composition develops at the outer border of the convective core. The atmosphere is convectively stable if the Ledoux criterion is met and unstable if the Schwarzschild criterion is met.^42,43 More explicitly, in order for convection to occur, the adiabatic temperature gradient should be smaller than the actual temperature gradient of the surrounding gas (Schwarzschild), which is given by the radiative temperature gradient if convection does not occur (Ledoux). We have examined the development of plumes as analogous to the overshooting observed in Schwarzschild-type instabilities.

Several people were involved in this work that is meant to honor Kenneth Showalter. I would like to acknowledge the work of Dr. Marcus Hauser of Otto-von-Guericke University in Magdeburg for calculating Marangoni factors and thermogravitational effects. Professor Bice Martincigh of the University of KwaZulu-Natal did the work shown in Fig. 5. Professor Tony Howes of the University of Queensland in St. Lucia worked on the data shown in Fig. 1. Special thanks and acknowledgements go to Matthew Eskew of Portland State University who did a lot of the experimental work that goes with Figs. 6–9. He also sketched Fig. 1. This work was supported by Research Grant No. CHE-1056366 from the National Science Foundation (NSF) and a Research Professor Grant from the University of KwaZulu-Natal, Cost Center No. P565.