Chemoconvection patterns in the methylene-blue – glucose system : Weakly nonlinear analysis

The oxidation of solutions of glucose with methylene-blue as a catalyst in basic media can induce hydrodynamic overturning instabilities, termed chemoconvection in recognition of their similarity to convective instabilities. The phenomenon is due to gluconic acid, the marginally dense product of the reaction, which gradually builds an unstable density profile. Experiments indicate that dominant pattern wavenumbers initially increase before gradually decreasing or can even oscillate for long times. Here, we perform a weakly nonlinear analysis for an established model of the system with simple kinetics, and show that the resulting amplitude equation is analogous to that obtained in convection with insulating walls. We show that the amplitude description predicts that dominant pattern wavenumbers should decrease in the long term, but does not reproduce the aforementioned increasing wavenumber behavior in the initial stages of pattern development. We hypothesize that this is due to horizontally homogeneous steady states not being attained before pattern onset. We show that the behavior can be explained using a combination of pseudo-steady-state linear and steady-state weakly nonlinear theories. The results obtained are in qualitative agreement with the analysis of experiments.


I. INTRODUCTION
When a solution of glucose contained in an open Petri dish and in the presence of methylene-blue is oxidized in a basic medium, the blue oxidized form of methylene-blue eventually becomes colorless (except in a thin layer at the free surface) [1].Gluconic acid is formed as a product and accumulates throughout the solution, particularly in the upper extent, where there is a continuous supply of oxygen from the atmosphere.This product forms a solution which is 0.044 g / cm 3 M more dense [2] than a corresponding equimolar solution of glucose, and the ensuing top-heavy distribution results in an overturning instability.Pons et al. [2,3] presented detailed experiments on the resulting patterns and time and length scales, and provided compelling evidence for the above instability mechanism.The reaction itself is well known [4][5][6] and its chemistry can be described using just two chemical equations [2] as a robust approximation to the full but complicated system [7].Bees et al. [8] formulated a minimal model of the hydrodynamic system and explored the linear stability of steady-state and pseudosteady-state (PSS) profiles.Linear theory predicts that the lowest value of the critical Rayleigh number for which the system first becomes unstable to a mode-1 solution occurs at a wavenumber approaching zero.This is not the more usual situation in the literature [9], although it is not exceptional.For instance, it occurs in convection with poorly conducting boundaries [10], for which convection develops in a layer of fluid heated from below such that the fluxes of heat at the top and bottom are constant and are not affected by convection.Also, a similar behavior can occur in bioconvection, which is ultimately driven by the biased motions of swimming microorganisms.Patterns reminiscent of convection, termed bioconvection, arise when motile microorganisms in suspension, with densities dissimilar to the fluid in which they swim, aggregate and consequently induce gravitational instabilities.For example, algae and bacteria may swim upward in response to gravity [11,12], light [13], or chemical gradients such as oxygen [14], and so accumulate at the top of the layer, or may accumulate in certain regions of the flow due to a balance of gravitational and viscous torques [12,15].Recent nonlinear analyses [16] of bacterial bioconvection about zero critical wavenumbers have highlighted the limitations of linear analyses of these systems.Furthermore, experiments [17] on bioconvection in suspensions of gyrotactic, swimming algae (for which linear theory predicts a competition between zero and nonzero critical wavenumbers [15]) reveal that the dominant pattern wavenumber varies significantly with time during the initial stages of the experiments in regions of parameter space close to the linear neutral stability curve.We shall show that analogous initial behavior in the chemoconvection system is not captured by the standard weakly nonlinear approach.
Amplitude equation based descriptions can usefully describe the dynamics of a system for a small, nonvanishing step above the linear neutral stability curve [9].In some cases the amplitude description is also a fair approximation of the full nonlinear system some distance from the bifurcation.Symmetry arguments allow for the construction of the form of the amplitude description for the system, although the explicit details must be obtained by other methods.For instabilities with a zero critical wavenumber, the precise amplitude description may, in general, be obtained by an appropriate expansion in terms of the long horizontal scale.Although this form of expansion is documented in the literature [10,16,[18][19][20][21] several nontrivial complications occur in the chemoconvective system, as we shall describe.
In Sec.II we describe essential experimental observations of the chemoconvection system.The amplitude description about a horizontally homogeneous steady state is derived in Sec.III, where the solution behavior is analyzed.Furthermore, we discuss the apparently conflicting predictions from PSS linear theory and the amplitude description.These results are reviewed in Sec.IV where we state some open problems.

II. EXPERIMENTAL OBSERVATIONS
An example of a typical chemoconvection pattern is presented in Fig. 1, observed from above [2].Figure 2 shows a plume of sinking fluid produced "artificially" from an oxygen source (a small hole in the lid of a completely filled cuvette; see Ref. [2] for details).This experiment, as well as other evidence [2], strongly supports the hypothesis of an overturning instability due to the formation of dense reaction products.
The essential reactions can be described by [2] 2MBH + O 2 → 2MB + + 2OH − , ͑1͒ where MB + is the blue oxidized form of methylene-blue, MBH is the colorless reduced form, GL is glucose, GLA is gluconic acid, and the rate of (1) is much larger than that of (2).
In Bees et al. [8] a minimal model of the system (fluid flow plus reaction-advection-diffusion) was constructed and the linear stability analysis of steady and pseudosteady profiles investigated.An important result of this is the prediction of the appearance at criticality of a pattern with wavenumber approaching zero.
In Pons et al. [3] the first quantitative analysis of the evolution of experimental patterns was documented.Therein, the dominant time and length scales were studied for sets of evolving patterns for which the experimental parameters were varied (particularly viscosity, pH, and layer depth).The overall behavior was similar in all cases and composed of three mostly identifiable stages.In the first, linear regime (very early stages of the evolution), the dominant wavenumber increases with time, which may possibly be explained with the PSS approximation (see later) and understood in terms of evolving neutral stability and linear growth curves.In the second (apparently transitionary) stage the growth of the dominant wavenumber slows and stops.In the third, and final, fully nonlinear regime a variety of behavior occurs (see Fig. 3).For instance, under certain conditions the dominant wavenumber monotonically decreases at long times or, in other cases, competition between modes is evident [3], resulting in self-sustained oscillations in the dominant wavenumber.The second and third stages cannot be satisfactorily explained using linear theory [8].
A transition from dot shaped patterns to line shaped patterns is observed when either the temperature or depth is increased.This transition may be a consequence of an exchange of stability between different configurations (dots and rolls) [22] or may be due to secondary physical mechanisms, such as surface tension effects or otherwise (note, however, that patterns also emerge and evolve in closed containers with oxygen permeable upper boundaries).FIG. 2. Image corresponding to a chemoconvective plume.The oxygen-rich plume sinks from a small hole in the lid of a completely filled cuvette of width 6.3 cm, height 6.7 cm, and thickness 1.3 cm.R is of order 1.The experiment is conducted at room temperature.
This paper aims to take a step in the direction of elucidating some of these nonlinear phenomena by means of an amplitude equation approach, although we shall restrict attention to two dimensions in this initial study.

III. DERIVATION OF THE AMPLITUDE EQUATION ABOUT A STEADY STATE
We begin from the following nondimensionalized set of governing equations [8] in two dimensions: where k is a unit vector directed vertically upward.Here, A, B, and ⍀ denote the concentrations of gluconic acid, the oxidized (blue) form of methylene-blue and oxygen, respectively, and u represents the fluid velocity.Furthermore, ␦ A and ␦ are diffusivity ratios relative to the diffusivity of the methylene-blue catalyst, and stand for reaction ratios, and R and S c correspond to the Rayleigh and Schmidt numbers, respectively.For completeness (see Bees et al. [8] for more details), we note that B and A have been nondimensionalized with W 0 , ⍀ with ⍀ ¯0, length with H ¯= ͱ D / k obs , and time with k obs −1 , where W 0 is the total (oxidized plus reduced) concentration of methylene-blue, ⍀ ¯0 is the concentration of oxygen at the upper surface, D is the diffusivity of MBH and MB + , and k obs is the effective rate constant [2] of reaction (2).Furthermore, it is worth noting that the Rayleigh number for this problem is defined as R = g⌬W 0 H ¯3 / d D, where g is the acceleration due to gravity, d is the viscosity, H ¯is the sublayer depth of reactants (as much of the instability is driven by reactions occurring in this surface sublayer), and ⌬ is the molar excess solution density of GLA, with respect to GL [2,8].
The first three equations model the reaction-advectiondiffusion of each reactant, and the last two equations are the Navier-Stokes equations for viscous flow with a gluconicacid-dependent density term, subject to the Boussinesq approximation, and an incompressibility condition, respectively.Solutions are required to exhibit zero vertical flux of A and B at the upper, z = 0, and lower, z =−d, boundaries such that Additionally, there is zero oxygen flux at the lower boundary,

͑9͒
and we impose fixed flux at the upper boundary,

͑10͒
where the constant C is chosen such that the concentration of nondimensional oxygen, ⍀, is equal to 1 at the surface before the onset of instability.It is clear that the imposition of this boundary condition constitutes a simplification of rather complex chemistry and transport mechanisms across a thin layer at the upper surface [23].
In this study we shall consider a two-dimensional system so that variations in the transverse y direction are zero.For a solid (no-slip) boundary, such as typically occurs at the lower boundary, u = ͑u,0,w͒ = 0, ͑11͒ while for a stress-free boundary, which typically occurs at the upper surface, The linear analysis predicts that the neutrally stable mode with the lowest value of the Rayleigh number (critical Rayleigh number R c ) occurs for a wavenumber of zero.However, a wavenumber of zero requires a system of infinite extent which cannot be realized experimentally.In a large finite system, with a Rayleigh number close to the critical one, one might expect to see a single plume (although the time taken for this instability to develop is long).Thus it is natural to investigate the behavior of the system with regard to a long horizontal length scale L. Therefore, we rescale the horizontal variable x with a small (but noninfinitesimal) parameter ⑀ = O͑L −1 ͒ Ӷ 1 [10], to give Furthermore, so that time evolution and leading order nonlinear effects (just) appear at the same order of ⑀ (in the (ϩ) corresponds to depth 0.68 cm and (ϫ) corresponds to depth 1.17 cm.The time shown is counted from the instant when the pattern appears in each experiment.Two experiments from many, as reported in [3].
resulting amplitude description), a systematic investigation of scales indicates that we need also to scale time and the Rayleigh number such that where R c = O͑1͒.Scaling otherwise leads to trivial solutions, linear analysis or systems of partial differential equations with complexity comparable to that of the original system.
The incompressibility condition (7) is automatically satisfied when u is written in terms of a stream function Additionally, we set = ⑀ ͑17͒ in order that the magnitude of the fluid velocity is at least one order higher in ⑀ than the rest of the variables; we impose zero flow before pattern onset.Note that the system necessarily incorporates at least two characteristic time scales as the mean concentration of gluconic acid always increases (in the time of any feasible experiment) subject to the availability of oxygen at the upper surface.Furthermore, note that although the horizontally homogeneous profile of gluconic acid never reaches a steady state, the variation of A about its spatial mean does.The vertical gradient, rather than the absolute value of A, drives the instability and so the equations are invariant to constant additions of A. Hence, the following ansatz generates a set of equations that takes account of the spatially averaged increase of A: where the superscript z indicates a vertical average and the subscript X indicates the partial derivative with respect to X. F͑X , T͒ is the amplitude function to be determined and represents a horizontally local excess or deficit of vertically averaged gluconic acid, and its presence in the definition of A is indicative of the aforementioned invariance.The functional dependence of on F X arises naturally, and is chosen ab initio for convenience, consistent with the governing equations.[Note that the variable 0 is not strictly part of the steady state (no-flow) solution as it has been scaled with ⑀, but labeled for convenience.] Hence, rescaling the equations, taking the curl of (6), substituting (18), and collecting terms of the same order in ⑀ gives the following systems of differential equations.At order ⑀ 0 The boundary conditions (8)-( 12) become Equivalently, for convenience at this order, we may fix the oxygen concentration directly at the upper surface by setting ⍀ 0 = 1, and maintain null flux at higher orders [see text after Eq. ( 10)].Equation ( 22) may be solved using these boundary conditions to give 0 = − 2z 4 − 3dz 3 + d 3 z 48 .͑25͒ For a realistic set of model parameters [8], = 248.0,= 21.9, ␦ A = 1.675, ␦ = 5.275, and d = 25, numerical solutions for A 0 , B 0 , and ⍀ 0 , and their derivatives, are presented in Fig. 4.

A. Solvability and solutions at order ⑀ 2
At second order, a solvability condition is required in order to remove secular (resonant) terms.This solvability condition can be obtained via Fredholm alternative theory or by inspection (due to uniqueness).
where M 0 ͑z͒ and N 2 ͑z͒ are defined in Appendix B. This equation is analogous to that obtained by Chapman and Proctor [10].
We thus obtain With this solvability condition and the realistic set of parameter values as above, solutions to Eqs. (B1)-(B8) may be computed (Fig. 5).However, A ˜2͑z͒ and A ˆ2͑z͒ are arbitrary up to the addition of a constant c, which consequently alters the profile of ˜2 and ˆ2 with the addition of c 0 ͑z͒ [see Eq. (B4)] and 2c 0 ͑z͒ [see Eq. (B8)], respectively.These factors can be absorbed into the function G͑X , T͒ which is determined at higher orders and so do not affect the amplitude equation in F. Hence, we are free to set A ˜2 z = A ˆ2 z = 0. Boundary conditions equivalent to these integral conditions are easier to implement numerically and can be obtained in the following manner.Equations (B4) and (B8) are integrated with respect to z to obtain relationships on the third derivatives of ˜2 and ˆ2 at the boundaries.For example, for a stress-free upper boundary and no-slip lower boundary we obtain Using these boundary conditions, solutions to the equations are presented in Fig. 5.

B. Solvability condition at order ⑀ 4
The solvability condition at order ⑀ 4 is determined in a similar manner to order ⑀ 2 .Integrating Eqs.(A1)ϩ(A2) ϩ͑ / ͒(A3) with respect to z over the fluid depth using the boundary conditions, we obtain, after some integration by parts, Here FIG. 4. Steady profiles at order ⑀ 0 (solid lines) and their derivatives (dashed lines).They are calculated using the following set of parameters: = 248.0,= 21.9, ␦ A = 1.675, ␦ = 5.275, S c = 2500.0(as standard but not used), and d = 25.The vertical step size is ⌬z = 25/ 1024.A 0 ͑z =−d͒ is fixed to zero for numerical convenience (we rename this field a 0 ), which has no further consequence in our calculations as only its derivative is used at higher orders.
Equation (39) is the required amplitude equation and models the evolution of the modulating function F. This amplitude equation has been obtained in other contexts (Rayleigh-Bénard convection between insulating boundaries [10] and Marangoni convection [19]).Equation (39) has the form of a conservation equation, expressing a symmetry in the system ( [9]; system invariant to homogeneous additions of A).
For the aforementioned realistic set of model parameters, we numerically calculate that R c = 4.893ϫ 10 −4 .Furthermore, we obtain the amplitude equation (39) with parameter values where is determined directly from the Rayleigh number of the system.Equation (39) represents the amplitude description for the evolution, in terms of the modulating function F, of finite amplitude perturbations about the horizontally homogeneous steady state.To obtain the amplitude equation in terms of the original space and time scales we substitute X = ⑀x, T = ⑀ 4 t and the definition of in the amplitude equation (39), such that

C. Steady solution of the amplitude equation
Chapman and Proctor [10] make a detailed study of the spatially periodic steady solutions of equation (39).Initially they consider k 3 = 0 and conclude that there exist infinitely many periodic solutions for an infinite system, of various spatial frequencies.These solutions can be obtained employing elliptic functions.Furthermore, qualitative changes in horizontally unconstrained solutions do not occur for different values of k 1 , or k 2 , as these can be absorbed into the variables and parameters with the transformations We choose the boundary conditions (i.e., such that there is zero horizontal flow and tangential stress at x =0, L).
For finite L, which quantizes the possible periodic solutions depending on the length of the box, there are n possible steady solutions differing in wavenumber (number of nodes) Finally, the main effect of including k 3 0 is to break the up-down symmetry of the solution.The larger the value of k 3 the more asymmetric will be the profile producing descending plumes [10].In our experiments up-down asymmetry is clearly evident: the down-welling centers of the dot or line patterns are always blue and the circulation cells' edges are always colorless.

D. Time-dependent solutions in a box
The time evolution of solutions of the amplitude equation can best be obtained numerically.For k 3 = 0, Chapman and FIG. 5. Steady profiles at order ⑀ 2 (solid lines) and their derivatives (dashed lines).They are calculated using the set of parameters shown in caption of Fig. 4.
Proctor investigated the stability of the steady-state solutions to small perturbations, both analytically and numerically (using a Lyapunov function representation).For k 3 0, stability results were also obtained for an appropriate initial value problem.The main conclusion from these results is that each mode is unstable to one of smaller wavenumber.This prediction is consistent with our experimental results for long times (e.g., see Fig. 3).This result is a genuine consequence of the nonlinear amplitude description and, therefore, cannot be explained using linear theory.However, the prediction is inconsistent with the experimental results in the early stages of the instability, which is likely to be due to the nonsteady nature of the horizontally homogeneous vertical profiles at onset.
In Figs. 6 and 7 we present the above behavior for a choice of parameter values.In particular, we can see that the wavenumber of the long-time state in the figure increases with ͑R − R c ͒ / R c .As before, to explain this one may remove the explicit dependence of the amplitude equation on ͑R − R c ͒ / R c by rescaling the length scale via Eq.(43).However, the scaled system length LЈ ϵ ͱ͓͑R−R c ͒k 1 / R c k 2 ͔ L then increases with ͑R − R c ͒ / R c .Since the disturbance is left to evolve to its natural state, more waves can fit into the larger rescaled box length.
In conclusion, the amplitude description is able to explain the qualitative experimental behavior at long times, but not at the initial stages.

E. PSS analysis versus the amplitude description
The PSS approximation [8] describes the linear stability of the system as horizontally homogeneous solutions evolve to steady state.It provides useful information as long as the system is close to the curve of (steady-state) neutral linear stability.For each time t, the PSS analysis predicts a neutral curve R͑k , t͒ and a linear growth rate curve ͑k , t͒ for modes of different wavenumber k, which slowly evolve to the steady-state curves.The neutral curve is seen to decrease ∀ k and the growth rate curve to increase.For a given set of parameters just above supercriticality at steady state, there will exist a time in the evolution of the horizontally homogeneous solutions where the system is critical (with a critical wavenumber of zero), with regard to the PSS approximation.Beyond this time the wavenumber of the mode with the greatest linear growth rate increases with time (the zero wavenumber always has zero linear growth rate [8]).Inherent in the PSS description is the assumption that any instability develops in a time frame that is quicker than that associated with the evolving profiles.There is experimental evidence to suggest that this is the case ([3]; profiles evolve over a time scale of an hour, and slower as steady state is approached, and patterns first emerge and then evolve in the linear regime over a time scale of 10 min).
However, as we observe above, the amplitude description indicates that each mode is unstable to modes of smaller wavenumber.
How can we reconcile these two views?Recall that the weakly nonlinear effects occur over a very long time scale, and we get a sense of this long term behavior from Fig. 3 (upper curve) and similar experiments [3].There are clearly time scale restrictions on the validity of each of the descrip- tions as discussed above which do not easily allow them to be merged.However, it seems reasonable to hypothesize that sufficiently close to criticality there is a transitionary period where PSS and weakly nonlinear effects act antagonistically.

IV. CONCLUSION
In this paper we have derived an amplitude equation for the chemoconvection system about a horizontally homogeneous steady state.If the system starts at steady state before instability onset, the amplitude description models the weakly nonlinear development and interaction of modes.The amplitude equation is of the form derived by [10] for which modes are unstable with respect to modes of smaller wavenumber (grain coarsening).
Previous linear studies [8] have highlighted the fact that a horizontally homogeneous steady state is never fully achieved before the onset of instability (which is also the situation with bioconvection, although never addressed theoretically).Linear analyses about these evolving profiles, termed pseudo-steady linear analyses, have demonstrated that the most unstable wavenumber increases with time from zero (at criticality).In this paper, we attempt to explain observations of experiments close to the curve of neutral linear stability using a blend of pseudo-steady-state linear analysis (leading to an increase in dominant wavenumber with time) and weakly nonlinear analysis about an initially steady, horizontally homogeneous solution (which predicts a long term reduction in the dominant wavenumber).In summary, at initial times the wavenumber grows, resulting from a slowly evolving vertical concentration profile that can be interpreted in terms of slowly evolving neutral stability and linear growth curves.In the long term, nonlinear terms control this linear growth and trigger a coarsening process.Clearly there must be a transitionary period in which the dominant wavenumber stops increasing and begins to decrease.Interestingly, experiments are suggestive of these three regimes.
An alternative interpretation might be valid for systems with Rayleigh number well above the critical one.For such experiments, a band of wavenumbers may grow rapidly associated with strong advection.This behavior could only be described by a fully nonlinear theory.However, it is possible that the states found in the weakly nonlinear theory are globally stable and that the transient solutions decay to them in the long term.
Future work should address the question of whether it is possible to construct a reduced description (such as an amplitude equation) that can incorporate this wealth of time scales, particularly in regimes where instability develops well before a steady state is approached, or whether the dynamics in each of these regimes are best separately described.
There is much experimental behavior in our system that cannot be currently explained, such as mode competition (e.g., self-sustained oscillations, observed in many experiments)-likely due to fully nonlinear effects well above the neutral linear stability curve.Future research should aim to explain these features.

FIG. 1 .
FIG. 1. Picture showing a chemoconvection pattern seen from above.The width of the image corresponds to 12.4 cm.The experimental conditions are as follows: temperature 19 °C, depth 0.59 cm, ͓NaOH͔ = 0.02 M, ͓methylene-blue͔ = 4.6ϫ 10 −5 M, ͓glucose͔ = 0.054 M, and R Ϸ 1.48.The image is taken 4800 s after pouring the fluid into the Petri dish.Darker regions correspond to the oxygen-rich fluid originally located on top of the layer.Lighter regions correspond to the fluid originally located below the blue subsurface layer.

FIG. 3 .
FIG. 3. Time evolution of the dominant wavenumber (k; periods per image width 12.4 cm) of chemoconvection patterns.Experimental conditions are identical to those in Fig. 1 except for depth.(ϩ)corresponds to depth 0.68 cm and (ϫ) corresponds to depth 1.17 cm.The time shown is counted from the instant when the pattern appears in each experiment.Two experiments from many, as reported in[3].
(32)-(35) into Eqs.(26)-(29) produces the system of equations given in Appendix B. Calculation of the adjoint operator of the homogeneous system, with the inner product of X and Y given by ͗X , Y͘ = ͐ −d 0 dz X • Y, is straightforward (Appendix A).The adjoint problem possesses the solution ͑A † , B † , ⍀ † ͒ = ͑1,1, / ͒.The solvability condition can be obtained by finding the inner product of the adjoint solution with Eqs.(26)-(28); directly integrating Eqs.(26)ϩ(27)ϩ͑ / ͒(28) with respect to z over the fluid depth gives an equation that determines the critical Rayleigh number for the instability:

FIG. 6 .
FIG. 6. Evolution in time of F͑x , t͒ from amplitude equation (42), with random initial conditions.The values of F͑x , t͒ have been normalized with respect to the final maximum and minimum values of F͑x , t͒, to 256 gray levels.Darker regions correspond to lower values of F͑x , t͒ and lighter regions to higher values of F͑x , t͒. x is represented horizontally and time t increases upward.The set of parameters in the amplitude equation for the model parameters used in Fig. 4 is as follows: k 1 = −1.67,k 2 = −93.85,k 3 = 2.67, k 4 = 0.285, ͑R − R c ͒ / R c =1, L = 1000, and total evolution time is 5000.Here, ⌬x = 1000/ 4095 and ⌬t = 0.2.