Nonlinear chemoconvection in the methylene-blue – glucose system : Two-dimensional shallow layers

Interfacial hydrodynamic instabilities arise in a range of chemical systems. One mechanism for instability is the occurrence of unstable density gradients due to the accumulation of reaction products. In this paper we conduct two-dimensional nonlinear numerical simulations for a member of this class of system: the methyleneblue–glucose reaction. The result of these reactions is the oxidation of glucose to a relatively, but marginally, dense product, gluconic acid, that accumulates at oxygen permeable interfaces, such as the surface open to the atmosphere. The reaction is catalyzed by methylene-blue. We show that simulations help to disassemble the mechanisms responsible for the onset of instability and evolution of patterns, and we demonstrate that some of the results are remarkably consistent with experiments. We probe the impact of the upper oxygen boundary condition, for fixed flux, fixed concentration, or mixed boundary conditions, and find significant qualitative differences in solution behavior; structures either attract or repel one another depending on the boundary condition imposed. We suggest that measurement of the form of the boundary condition is possible via observation of oxygen penetration, and improved product yields may be obtained via proper control of boundary conditions in an engineering setting. We also investigate the dependence on parameters such as the Rayleigh number and depth. Finally, we find that pseudo-steady linear and weakly nonlinear techniques described elsewhere are useful tools for predicting the behavior of instabilities beyond their formal range of validity, as good agreement is obtained with the simulations.


I. INTRODUCTION
Hydrodynamic instabilities arising from chemical reactions ͓1-4͔ can have a fundamental impact on the overall reaction rates of many processes ͓5͔ and are compelling examples of self-organized pattern formation.A deep understanding of the role of pattern formation could benefit a wide range of industrial applications, and may lead to optimal methods of production and control.Yet in general, we do not yet have a complete picture of the pattern forming process and can say much less about nonlinear evolution of the pattern and its potential impact on aggregate reaction rates.
In contrast, nonhydrodynamic pattern formation in chemical and biological systems has been studied in great theoretical detail in recent decades.The most well known are the Belousov-Zhabotinsky class of reactions, possessing targets and scroll waves, and autocatalytic systems exhibiting Turing patterns ͑see ͓6͔, and references therein͒.On the other hand, the well-known Rayleigh-Taylor and Rayleigh-Bérnard problems, and many variants thereof, have also been well documented and theoretically explored ͑e.g., ͓7,8͔͒, as well as general hydrodynamic instabilities too numerous to mention ͓9͔.As many chemical engineering processes involve reactions in fluids and products with a variety of physical properties it would appear that detailed theoretical studies of hydrodynamic instabilities caused by reactions are long overdue.
Biologically induced hydrodynamical instabilities ͓10,11͔, however, have been studied in theoretical detail ͓10,12-24͔.In these systems, the swimming behavior of single-celled micro-organisms, in response to light, chemical gradients, or viscous and external torques ͑such as gravity͒, leads to aggregations of typically relatively dense cells from which an instability can result.In such biological systems the "energy" of the instability is provided by the swimming behavior of the individual cells, whereas energy typically is provided at the boundaries in classical convection problems ͓7͔.Recently, interest in mechanisms that induce hydrodynamic patterns in chemical solutions, particularly at interfaces ͓2-4,15-19͔, has increased.Theoretical modeling and linear and weakly nonlinear stability analyzes ͓20,21͔ have been reported but a mechanistic understanding of full nonlinear behavior is lacking.In this paper, we explore comprehensive two-dimensional nonlinear simulations.Due to the complex nature of general chemical-hydrodynamic systems, we choose to study an underlying chemical system that is experimentally and theoretically simple; we investigate the well-known methylene-blue-glucose system, also called the "blue-bottle experiment" ͓22͔.The complete chemistry of this system is actually fairly complex, but studies have shown it to be amenable to significant simplification ͓4͔.
The main aim of this study is to directly compare twodimensional numerical solutions of a recent model of the methylene-blue-glucose system with experimental observations of chemoconvection patterns and previous analyses.To achieve this aim it is necessary to explore the sensitivity of the system to the form of one crucial boundary condition, how oxygen enters the system; the upper boundary condition presented analytical issues in previous analyzes and it affects the solution behavior qualitatively.Furthermore, we shall explore parameter space and the aspect ratio of the container concentrating on temporal and spatial scales.
In the next section we describe preliminaries of the methylene-blue-glucose system and summarize experimental and theoretical results from previous work on pattern formation in this system.In Sec.III we describe the simplified spatial model, and in Sec.IV the numerical scheme.In Sec.V we present the results.We discuss comparisons with experiments, and linear and weakly nonlinear theories in Sec.VI and go on to calculate the impact of pattern formation on aggregate rates of reaction.Finally, we present conclusions in Sec.VII.

II. THE METHYLENE-BLUE-GLUCOSE SYSTEM
When a mixture of glucose ͑GL͒, sodium hydroxide ͑NaOH͒, methylene-blue and water is shaken with oxygen ͑O 2 ͒ and poured into a Petri dish that is open to the atmosphere a set of chemical reactions ensues ͓17,23͔ ultimately producing gluconic acid ͑GLA͒.The essential reactions can be represented robustly ͓4͔ by 2MBH + O 2 → 2MB + + 2OH − , ͑1͒ GL + MB + + OH − → MBH + GLA, ͑2͒ where MBH and MB + are the colorless ͑reduced͒ and the blue ͑oxidized͒ forms of methylene-blue, respectively.The rate of oxidation, k 1 , in Eq. ͑1͒ is much greater than the rate of reduction, k obs , in Eq. ͑2͒ ͑measured by Pons et al. ͓4͔͒.
Thus in the presence of oxygen the mixture turns blue as MBH is quickly oxidized throughout the fluid ͑k 1 ӷ k obs ͒.However, once in the container and at rest, the mixture slowly returns to the colorless form as oxygen is used up, except at the upper open boundary where there is a supply of oxygen.Here, both reactions ͑1͒ and ͑2͒ proceed unheeded and the reactants diffuse, generating a positive vertical gradient in GLA concentration.Pons et al. ͓4͔ established that this product forms a solution that is 0.044 g / cm 3 M more dense than an equimolar solution of glucose.Therefore, a gravitationally unstable layer is formed at the upper surface that can lead to an overturning instability.Tens of minutes after the solution is poured into the Petri dish, a pattern of blue sinking fluid in a colorless background typically emerges and evolves over the course of many hours ͑see examples of these patterns in Fig. 1, as viewed from above͒.
Pons et al. ͓4͔ presented conclusive evidence for the nature of the instability mechanism, supported by two simple experiments.In the first, they found that similar patterns can be produced in a container with an oxygen permeable lid ͓Fig.1͑c͔͒, implying that surface tension effects are not necessary to produce the instability.In the second, a small hole was made in the lid of an oxygen impermeable container resulting in a thin, blue, descending plume of oxygen rich fluid, indicating that density variations are sufficient to drive the instability.They went on to observe the qualitative dependence of the instability on pH ͑affecting reaction rates͒, glucose, catalyst concentration, fluid depth, and temperature.
The first comprehensive quantitative analysis of the pattern growth and wavelength was described in Pons et al. ͓15͔ subject to a variety of experimental conditions.In particular, they explored pattern length and time scales due to variations in viscosity, pH, and fluid depth.After the unstable vertical density profile evolves to such an extent that the initial instability arises and a pattern is formed, the dominant wavenumber is seen to increase quickly.This rapid growth in dominant wavenumber is later retarded: the dominant wavenumber either is maintained, slowly decreases, or oscillates ͑with a fairly well-defined period͒ for long times, suggesting a regular sequence of mode interactions.
Bees et al. ͓20͔ developed a reaction-advection-diffusion system to model the pattern formation mechanism.Linear stability analysis was conducted from a horizontally homogeneous, no-flow, steady-state concentration profile with a fixed oxygen concentration upper boundary condition ͑hereafter referred to as FCBC͒, using independently measured parameter values ͓4͔.They found that the critical wavenumber is zero due to symmetries of the governing equations ͓21͔.͑The critical wavenumber is associated with the smallest Rayleigh number, R c , on the neutral linear stability curve.For Rayleigh numbers just greater than R c , a band of unstable modes exist with small wavenumber.͒Nevertheless, the fastest growing linear mode for a Rayleigh number above R c has a finite wavenumber; this wavenumber and its growth rate have been favorably compared with experimental results ͓15,20͔.
However, simple linear stability analysis based on steady profiles was not able to explain much beyond simple attributes of the initial pattern.For instance, it could not explain the onset time at which the pattern emerges or elucidate the initial quick increase of the dominant wavenumber after the instability arises.Moreover, it is clear from experiments that a no-flow, steady-state concentration profile is generally far from attained before pattern onset.Bees et al. ͓20͔ introduced a pseudo-steady state approximation to address these aspects.As the vertical chemical concentration profiles asymptote to the steady state, linear stability analyzes were conducted about slowly evolving concentration profiles.Even though the instability may be on the verge of being induced ͑with vanishingly small growth rates͒, the vertical profile will continue to evolve and become more unstable.This neutral curve can be considered to be a slow function of time and to move through a line of constant Rayleigh number ͑given by the physical properties of the reactive fluid͒, presenting a wider band of unstable modes with an increasing value of the most unstable wavenumber as the density profile steepens.The validity of the pseudo-steady state linear stability analysis is of course restricted by several constraints to a small range of parameters in the model.However, results appear to compare remarkably well with experiments, both qualitatively and quantitatively, outwith their expected range of validity ͓20͔.
In Pons et al. ͓21͔, a weakly nonlinear extension of the model was presented for a fixed oxygen flux upper boundary condition ͑hereafter referred to as FFBC͒ to establish the behavior of the system beyond the initial instability.The linear analysis closely follows that for the FCBC, and the resulting nonlinear amplitude equation is analogous to that for convection with insulating walls.Weakly nonlinear analysis indicates that dominant wavenumbers should decrease in the long term, thus presenting a trend opposite that seen for pseudo-steady linear stability analysis with a FCBC.However, such analysis suffers from similar problems as the steady linear analysis as steady state profiles are never attained before the instability occurs.In Pons et al. ͓21͔ it was suggested that a combination of pseudo-steady linear stabil-ity analysis and weakly nonlinear analysis from a steady state could reproduce much of the early behavior observed in experiments.As well as comparing our numerical solutions with experimental observations in later sections, we shall also discuss the range of validity of the above theoretical descriptions in an attempt to reconcile these seemingly different views.

III. MODEL EQUATIONS
Throughout this work, we shall consider two-dimensional solutions, such that solutions do not depend on y and the fluid velocity u = ͑u ,0,w͒.We shall begin with the nondimensionalized model as presented in Bees et al. ͓20͔ which consists of the following set of governing equations: where k is a unit vector pointing upwards.A, B, and ⍀ denote the nondimensional concentrations of gluconic acid, oxidized ͑blue͒ form of methylene-blue, and oxygen, respectively.␦ A and ␦ are diffusivity ratios and and are reaction ratios.R and S c represent the hydrodynamical nondimensional parameters and correspond to the Rayleigh and Schmidt numbers, respectively.Length has been nondimensionalized with H ¯= ͱ D / k obs , where D is the diffusivity of MBH and MB + , and time with t d = k obs −1 ͑recall that k obs represents the rate of the slow reaction͒.The chemical Rayleigh number is defined as R = g⌬W 0 H ¯3 / d D, where W 0 is the total ͑oxidized plus reduced͒ concentration of methyleneblue, g is the acceleration due to gravity, d is the dynamic viscosity, and ⌬ is the molar excess solution density of GLA, with respect to GL ͓4,20͔.Typical values of these quantities are noted in Table I.
The first three equations represent the reaction, advection, and diffusion of the three main chemicals: GLA, MB + ͑MBH can be determined from the initial homogeneous concentration of methylene-blue minus MB + ͒, and O 2 .The fourth equation is the Navier-Stokes equation with a negative buoyancy term due to GLA ͑Boussinesq approximation͒.And the fifth equation is the incompressibility condition.
We shall impose no-slip and no-flux boundary conditions at x =0,L: and u = ͑u,v,w͒ = 0 at x = 0,L.͑8͒ For a typical numerical experiment with no-slip at the lower boundary and stress-free at the upper boundary we impose, respectively, ͑u,w͒ = ͑0,0͒, and w = ‫ץ‬u ‫ץ‬z = 0. ͑9͒ Vertical no-flux boundary conditions are applied for A and B at the upper and lower surfaces such that

͑10͒
and at the lower boundary only for oxygen flux,

͑11͒
In Bees et al. ͓20͔ the concentration of oxygen at the upper surface was fixed to a value of 1 ͑nondimensional units; FCBC͒.But, in Pons et al. ͓21͔ a fixed oxygen flux upper boundary condition ͑FFBC͒ was investigated.Both boundary conditions are, at best, approximations of natural oxygen transport at the surface.The FCBC imposes a large flux if oxygen poor fluid is advected close to the surface, whereas a FFBC would clearly be inappropriate in the absence of an oxygen sink ͑also ͓17͔, presents evidence that the oxygen flux depends on the stirring rate of the system͒.
In this paper, we shall explore more general boundary conditions at the upper surface than considered elsewhere ͓20,21͔.We use mixed boundary conditions where the constants ␤ c and ␤ f are chosen such that we can obtain fixed concentration ͑␤ f =0͒, fixed flux ͑␤ c =0͒, or mixed type boundary conditions at the upper surface.For a FFBC, we choose ␤ f such that the concentration of nondimensional oxygen, ⍀, is equal to 1 at the surface before the onset of instability.We may rewrite condition ͑12͒ in the form of a naïve model of flux across a thin mixing layer, such that where ⍀ sat =1/ ␤ c is the oxygen saturation concentration and ⌬ = ␤ f / ␤ c is the thickness of a thin absorption zone ͑nondimensional units͒.

IV. NUMERICAL TECHNIQUES
Here we describe the numerical scheme to simulate solutions of the full nonlinear model, the standard independently measured parameter values used in the simulations, and methods to calculate statistics and Fourier spectra from the numerical results.

A. Numerical methods
We solve the evolution equations using a time-splitting method with an improved boundary condition for the pressure as described in ͓24͔.The time integration scheme is second order accurate and is based on a modified Adams-Bashforth formula ͓24͔; we treat the linear diffusion terms implicitly and the nonlinear terms explicitly.For the spatial discretization we use a Chebyshev collocation pseudospectral method ͓25͔.The advantage of using spectral methods for the spatial discretization is the high accuracy that can be achieved if the solution is sufficiently smooth; this leads to very small numerical diffusion that precludes the propagation of errors during the time evolution.A similar code has been used successfully for the simulation of binary fluid convection in large aspect ratio containers ͓26͔ showing excellent agreement with experiments.In all cases the time step and the number of collocation points used was adjusted until the solutions converged.Typically, we used 480 collocation points in the x direction and 48 collocation points in the z

B. Parameter values
The standard nondimensional parameter values are chosen to represent a perceived normal experimental setup, and are tabulated in Table II.For fixed-flux boundary conditions we choose the value of the flux such that the concentration of oxygen is equal to 1 at the surface before the onset of instability, in order to make direct comparisons with the fixed concentration boundary condition.This leads to a boundary concentration coefficient ␤ c of zero and a boundary flux coefficient ␤ f of 5.48, giving a flux of 0.182.We later explore the impact of varying ␤ c and ␤ f , before focusing upon the mixed condition values ␤ c = 1 and ␤ f =3.
In all cases, we employ the same initial conditions at t =0: A =0, B =1, ⍀ =1, ͉u͉ = 0; initially, the medium is quiescent, homogeneous, and blue, and it is saturated with oxygen.

C. Statistics and Fourier transforms
Taking into account the large quantity of data produced by each numerical experiment we define several global quantities which can help to simply describe the solution behavior.In particular, we make statistics of the concentration distribution of each compound and of the modulus of velocity, in the vertical and horizontal directions and for the whole system.These will be defined for the oxygen concentration below and it is clear how the definition is extended to other variables.
The oxygen average for the whole system is defined as computed using Clenshaw-Curtis weights.The vertical and horizontal averages, ͗•͘ z and ͗•͘ x , are calculated, similarly, restricting the integrals to the vertical and the horizontal directions.͑These averages are distributed in space in the horizontal and vertical directions, respectively.͒Also, we define the oxygen variance statistic, Var͑⍀͒, as

͑15͒
which is numerically computed in a similar manner to ͗⍀͘ above.Var͑⍀͒ gives a quantitative measure of the spatial inhomogeneity of the field.Similarly, we can define the vertical and horizontal analogs, Var z ͑⍀͒ and Var x ͑⍀͒.
͗B͘ z is what one typically observes in experiments, i.e., the depth-averaged pattern of methylene-blue from above.Therefore we study the evolution of the dominant wavenumber by Fourier transforming ͗B͘ z : for each time step interpolate ͗B͘ z ͑x͒ at equispaced values in x, apply the Hann windowing function and employ the one-dimensional FFT algorithm ͑e.g., see ͓11͔͒.The dominant wavenumber k max is defined as that associated with the largest value in the power spectrum.

A. Solutions for the standard parameter values with a fixed flux boundary condition
We begin by presenting snapshots, Figs.2͑a͒-2͑c͒, of the three concentration fields with the fluid velocity superimposed for a typical simulation output with the standard parameter values and FFBC ͑i.e., ␤ c = 0 and ␤ f = 5.48͒.This boundary condition allows for analytical results that can be compared directly with the simulations.We also present plots of Var x ͑A͒ and Var x ͑B͒ in Fig. 2͑d͒ to illustrate the form of the inhomogeneous part of the vertical concentration fields.In the very early stages, these can be favorably compared with plots of the perturbation fields from weakly nonlinear analysis ͑see Fig. 5 in ͓21͔͒.The long time profiles of Var x ͑A͒ and Var x ͑B͒ are also presented to illustrate that variations in B are towards the bottom of the fluid layer whereas variations in A extend to upper layers before decreasing close to the free surface.In Fig. 3, graphs are presented for the evolution in time of the average concentrations of gluconic acid ͗A͘, methylene-blue ͗B͘, oxygen ͗⍀͘, and modulus of velocity ͉͗v͉͘, where the inset of each graph depicts the initial stages of the evolution.Furthermore, we plot vertically averaged quantities ͗A͘ z , ͗B͘ z , ͗⍀͘ z , and ͉͗v͉͘ z in Fig. 4. Initially the mixture is homogeneous and so methyleneblue is oxidized at the same rate throughout the layer of fluid.However, for the FFBC by definition there is a given influx of oxygen at the upper surface from the atmosphere, leading to a thick oxygen depleted layer at the bottom.Hence methylene-blue is preferentially and rapidly oxidized at the upper surface and slowly reduced generating the dense product gluconic acid A. At a time of approximately 18 units ͑18 t d =18/ 0.0042Ϯ 240 s Ϸ 71Ϯ 4 min͒ the evolving balance of accumulation and diffusion of A at the upper surface leads to an overturning instability.Advection ensues and highly inhomogeneous oxygen concentrations are generated, as can be observed by a sharp increase in Var͑⍀͒ ͑inset Fig. 3͒.However, this oxygen distribution relaxes to a generally thin subsurface layer of oxygen that is mostly stable to advective effects but with a thin plume of entrainment in regions of downwelling fluid ͑see vertically averaged quantities in Fig. 4͒.The behavior of the oxidized methylene-blue B mostly follows that of oxygen, except that it is continually advected downwards from the subsurface layer, interacting with the lower boundary and generating the striking blue dots and rolls that can be observed in experiments for long times.Initially, B forms a thick band ͑thicker than oxygen͒ at the upper surface.After the overturning instability arises, thick plumes of B descend to the lower boundary, mixing the fluid and being advected back upwards, almost resulting in a homogeneous distribution of B. Excess B is quickly established once more at the upper surface and advected downwards to create a new array of blue plumes.There are clear oscillations in this behavior before more stable arrays of plumes are established.Figure 4 demonstrates this early oscillation in the amplitude of the pattern.
The difference between the distribution of B and oxygen is due to the relatively slow reaction rate of the reduction process.It is the inhomogeneity of the gluconic acid distribution, however, which drives advection; it can be observed in Fig. 2 how vertical gradients in A quickly arise and lead to thin plumes transporting excess A at the surface to broad accumulations at the bottom.
After transients have diminished, Fig. 3 appears to present an inherently stable picture of the spatially averaged quantities.Long term behavior of the averaged quantities indicates that steady states are approached for ͗B͘, and that ͗A͘ in- creases linearly with time.In fact, we can average Eq.͑3͒ obtaining ‫ץ‬ t ͗A͘ = ͗B͘, ͑16͒ and Eqs.͑4͒ and ͑5͒ to obtain and where ͗•͘ x means horizontal average.Equation ͑16͒ provides the rate of growth of ͗A͘ if ͗B͘ is known.Then, considering a steady state in the averages ͗B͘ and ͗⍀͘ gives

͑19͒
Setting a FFBC provides a steady state value for ͗B͘ of ͗B͘ = ͑␦/d͒‫ץ‬ z ⍀͑0͒ = 0.452, ͑20͒ which agrees ͑error generally less than 1%͒ with Fig. 3, and provides a ͑poor͒ positive lower bound for The fact that ͗B͘, and so the rate of increase of ͗A͘, numerically reaches a steady state value as directly calculated above is quite remarkable.In fact it is the same steady state value independent of whether convection arises or not.The evolution with time of ͗B͘ varies for small time intervals as the assumed steady state for ͗⍀͘ occasionally is broken, but resumes its convergent approach as steady state conditions for ͗⍀͘ are reattained ͑note that the "temporary" steady state value of ͗⍀͘ is not reattained; readers should examine closely Fig. 3͒.Such shifts in ͗⍀͘ are due to plume interaction events.Descending plumes in the initial stages of the instability advect oxygen down to oxygen-poor ͑colorlessmethylene-blue͒ regions.This oxygen is quickly used at the boundaries of the plumes before a balance between advection and use of oxygen is attained in the regions surrounding the plumes.
An implication of Eq. ͑20͒ is that the solution will not approach a temporary mean steady state for ͗⍀͘ if the parameters are chosen such that ‫ץ␦‬ z ⍀͑0͒ Ͼ d, as ͗B͘ is bounded above by 1.As an aside, setting the value of the fixed flux to three times previous values leads to a prediction that ͗B͘ = 1.35, which is clearly not possible.However, the simulations indicate that ͗B͘ approaches the value 1 after circulation is established.Therefore we re-examine our assumptions to find that a steady state in ͗⍀͘ will not be reached.Thus we find that with ‫ץ‬ z ⍀͑0͒ = 0.546, which is indeed observed in the simulations ͑error less than 1%͒.With this high value of oxygen flux, solutions are generated with a thick band of GLA at the upper surface and narrow band of less dense mixture at the bottom, leading to thin ascending plumes of GLA deficient fluid ͑not shown͒.
Let us now return to the example with ‫ץ‬ z ⍀͑0͒ = 0.182.The apparent steady states discussed above hide the variation in the solutions of the full system seen in simulations.In particular, one can observe a wealth of plume interactions, which in this FFBC case result in cell annihilation events and a final single circulation pattern occupying the whole of the container leading to a large horizontal gradient in ͓GLA͔.Fig. 3 demonstrates that although ͗B͘ tends to a steady state, and ͗A͘ increases as in Eq. ͑16͒, ͗⍀͘ and ͉͗u͉͘ do not reach steady states: ͗⍀͘ and ͉͗u͉͘ adopt new larger values after each plume collision ͑see Fig. 4 for two-plume merging events in the early stages of the pattern evolution͒.
We can observe how the wavenumber of the dominant mode of B, k max , varies with time in Fig. 5͑a͒.In Fig. 5͑c͒ we also plot the power spectrum of B from which we extract the maximum amplitude, as in Fig. 5͑b͒, and plot the normalized power spectrum with time in Fig. 5͑c͒.From these figures, and other simulation output, it is clear that there is a large increase in k max in the very early stages of the instability, followed by a gradual decline with some oscillations.In Fig. 6 we see that plume interaction events are direct and result in fewer plumes and thus a sudden decrease in wavenumber.͑We shall see later that this is not the case for the FCBC.͒We shall focus on two significant events at times Ϸ130 and 710 which are just discernible in Fig. 3 ͑curve 2͒ as steps up in ͗⍀͘ and ͉͗v͉͘, and blips in ͗B͘.These events can be matched with behavior in Fig. 4. The first event at time 130 consists of a merging of two large plumes to form a thick plume of downwelling blue fluid, whereas the second event at time 710 represents the collision of the thick central plume with a boundary.One may also observe in Fig. 4 small amplitude, high frequency wiggles on the plumes.In movies of the simulation results it is clear that these are not numerical artifacts but small traveling undulations which appear close to wide plume structures.
We note here a limitation of the simplified model kinetics that we have employed.It does not make sense to produce more ͓GLA͔ than the original ͓GLU͔.Essentially, the approximation assumes that ͓GLU͔ is constant when in fact it should diminish in response to increasing ͓GLA͔.We can estimate the time at which ͓GLA͔ = ͓GLU͔ and the model is strictly invalid.For times approaching this time, the results should be treated with caution.Initially, ͓GLU͔ is set to 0.054 M. Expressed in nondimensional units this is 0.054/ 4.6ϫ 10 −5 = 1174.The simulations do not pass this value, but ͗A͘ does approach half this value over a nondimensional time of 1000, corresponding to 1000/ k obs s = 66 h.However, experimental results do not exceed 5 h in ͓15͔, so the simulation times are well beyond what are required, although it makes sense to investigate simulation extremities.The influence of the reduction in ͓GLU͔ throughout the course of a normal experiment or simulation may modify the solutions quantitatively in the long term.
In what follows, we systematically explore the parameter space around the standard parameter values, and the boundary conditions.To illustrate the full range of solution behavior, we concentrate on representative changes of ͑i͒ upper oxygen boundary conditions, ͑ii͒ Rayleigh number, R, and ͑iii͒ depth to sublayer depth ratio, d.We also briefly tabulate dependence on other parameters such as ratios of reaction rates and diffusivities in Table III.

Fixed concentration
First we consider the case for a FCBC, which can be achieved by setting ␤ c = 1 and ␤ f = 0.At first glance, this would appear to be the most natural approximation to the real boundary condition for oxygen at the surface, but implies that oxygen present in the air takes no time at all to cross the gas-fluid interface, to the specified concentration.Moreover, it may be viewed as a much harsher restraint than the FFBC as it states that the concentration of oxygen at the free surface is not perturbed by advective or reactive effects, which can lead to very high spatial oxygen gradients, perhaps unphysically.In contrast, the FFBC, although directly prescribing a constant supply of oxygen, allows for oxygen concentrations to vary at the upper boundary and does not lead to such high gradients.It is apparent, however, that there are modeling issues with both types of boundary condition.
For FCBC it is immediately clear from the simulation results that oxygen is not at all entrained in the plumes, de- FIG. 5. ͑Color online͒ Variation of the dominant wavenumbers with time.͑a͒ the wavenumber with the greatest growth rate, k max , ͑b͒ the maximum amplitude of the power spectrum, ͑c͒ power spectrum ͑left͒ and power spectrum normalized with respect to its maximum at each time ͑right͒.The normalized power spectrum is represented only for times at which its maximum is larger than 10 −5 .͑a͒ and ͑b͒ are for standard parameter values with boundary conditions ͑1-5͒ as for Fig. 3, and ͑c͒ is for standard parameter values with a fixed flux oxygen boundary condition ͑FFBC͒.parting from the behavior for the FFBC.We plot in Figs.7͑a͒ and 7͑b͒ the early and long time behavior of ͗A͘ z − ͗A͘ and ͗B͘ z and, in Fig. 7͑c͒, a snapshot of A − ͗A͘ and B concentrations.The behavior is clearly different to that observed for the FFBC, with plumes tending to separate rather than merge.Common to both mechanisms is that wavelengths increase; for the FCBC, plumes tend to drift towards the edge of the container and annihilate at these boundaries.The plumes also tend to exhibit large amplitude wiggles with a period of approximately 150-200 time units.There are also a few large plume-translation events, such as at time Ϸ750 which can be associated with minor blips in ͗B͘ but better by a larger step and a blip in ͉͗v͉͘ in Fig. 3, respectively.Pairs of closely sited plume-antiplume "doublets" can be seen to rapidly translate and appear to cross the path of other doublets ͓see Figs.7͑a͒ and 7͑b͔͒.This is distinctly different behavior to the merging events that occur for the FFBC, and it seems remarkable that such a variety of behavior can be observed for such a simple change in boundary condition.There are no major steps in ͗⍀͘ in contrast to the FFBC.It is also apparent in the snapshots of Fig. 7͑c͒ that circulation rolls have greater up-down symmetry than the FFBC, with blue downwelling MB + rich plumes of a similar width to clear upwelling MBH rich zones.However, the solutions lack leftright symmetry, with the asymmetry directing the drift of meandering plume-antiplume doublets, typically outwards to the boundaries.

Mixed boundary condition
This difference in behavior between FFBC and FCBC, with the physical problems associated with each, begs the question as to whether an intermediate scenario may shed some light.Thus we consider a mixed boundary condition.Typically, solutions are obtained with elements of both fixed flux and concentration.By setting ␤ c = 1 and ␤ f = 3 a regime similar to the case of FFBC is obtained, with thick downwelling MB + rich plumes that merge, except that the plumes merge faster and almost simultaneously, and the highfrequency wiggles previously seen on the plumes have greater amplitude and larger period ͑roughly 12 time units͒.Furthermore, observation of the simulation results reveals interesting plume development and structure as can be observed in Fig. 8, and a treadmill of smaller undulations that get translated and enveloped by larger plumes before they have a chance to develop separate circulation cells.Oxygen is entrained in the plume structures.A state is quickly ob-tained in which the large central plume migrates to a boundary leaving a central upwelling region with downwelling at the boundaries.
Increasing ␤ f to 9 provides further new behavior.This region of parameter space prescribes a constant total contribution of oxygen concentration and a small flux.Alternatively, it may be viewed as physically modeling a broad absorption zone ͑⌬ = ␤ f / ␤ c = 9 nondimensional units͒.In the early stages, well-formed, thin downwelling plumes merge.At later stages and for long times two rather unstable-looking plumes gobble up smaller plumes in a state of nervous coexistence, drifting around in their allowed space but never desiring to get sufficiently close to merge.This leads to a large horizontal gradient in GLA, with higher concentrations in the well-mixed central zone housing the coexisting plumes.Of note is that oxygen is only very weakly entrained in the plumes and that large ͑but not all͒ plumes repel one another, closer in behavior to the case of FCBC.
Reducing ␤ f to 1 or 0.2 provides at first a regular array of plumes that quickly translate to the right or left so that they preferentially fill one half of the horizontal domain, eventually merging and translating to the boundary to form one large circulation roll, generating fairly constant spatial gradients in ͗•͘ z quantities.The behavior in these mixed boundary condition cases is generally significantly different to the FCBC described above, and suggests that the FFBC better reflects the behavior of more complicated boundary conditions.In the absence of reliable measurements on realistic values of ␤ c and ␤ f , we assume a mixed boundary condition with ␤ c = 1 and ␤ f = 3 in what follows.

C. Varying the Rayleigh number
Here, we describe the influence of the Rayleigh number on the evolution of the system.This parameter is important in the linear stability analysis of the system.Consider a reduction of the Rayleigh number from the standard parameter values.In general we find that the magnitude of any circulation that may arise decreases with R. Setting R =10 −1 we find that an instability still occurs with a similar wavenumber to the standard parameter values.However, it can be seen in Fig. 9 that the plumes translate much less rapidly and merge at times much later than for the standard parameter values.Furthermore, the plume merging events are less smooth: the plumes gradually drift towards each other followed by a very rapid coalescence once they are within a certain minimum distance.The solutions head towards the final solution for the standard parameter values, with the wavenumber decreasing monotonically.Oxygen is still entrained, but less so than for the standard parameter values.Of particular note is that there appears to be only a narrow range of allowable plume separation distances, perhaps associated with the smaller range of growing modes from linear theory.Further reducing R to 10 −2 , we find very smooth curves for the ͗•͘ z quantities as a function of time, with very few of the oscillations, steps, or blips as seen before.All quantities converge rapidly to steady states, consistent with the observation that no plume merging or translation events occur.In fact the long time behavior of the solution is a stable, sinusoidal array of plumes with a wavenumber commensurate with that of the initial instability.This leads to the suggestion that the results display behavior very similar to that predicted by linear theory, and that this value of the Rayleigh number should be close to the neutral stability curve ͑see Sec.VI A for a comparison with predictions from linear theory͒.A wavenumber of 5 in an x region of 240 gives k = 0.021 ͑2 s.f.͒ per unit, or k = 0.13 per 2, consistent with the linear theory prediction ͓20,21͔ of approximately k = 0.1 for R =2R c ͓which gives the scaled Rayleigh number r = ͑R − R c ͒ / R c =1͔, both for fixed concentration ͑R c = 5.2ϫ 10 −4 ͒ or fixed flux boundary conditions ͑R c = 4.8ϫ 10 −4 ͒, rather than the mixed boundary conditions, albeit for slightly different parameters.
Finally, setting R Յ 5 ϫ 10 −4 does not lead to instability, consistent with linear theory ͓20,21͔.Furthermore, the simulations indicate that the amplitude of nonlinear solutions decays with proximity of R to R c in a fashion consistent with a supercritical bifurcation.

D. Dependence on depth
The stability of plumes in deep and narrow chambers is of great interest, but beyond the scope of this paper ͑see discus- sion͒.On the other hand, if d Ӷ L, then we find that increasing L has minimal effect, so we consider varying depth d with L fixed.
Decreasing depth from the standard parameter values by setting d =12 ͑24 collocation points in the z direction͒ leads to results that exhibit a monotonic decrease in k.The initial instability generates circulation cells with small wavelengths comparable with the depth d, that transport ⍀, MB + , and GLA very rapidly throughout the fluid.Thus the initial wavenumber is larger for this shallow system, but is reduced very quickly.The circulation cells merge rapidly with their neighbors.At the end of the simulation time of 1000, the system ends up with a single centrally located plume that locally transports fluid upwards, with downwelling currents at the edges.
If we increase depth from the standard parameter values by setting d = 35 we find, as one might expect by analogy with the analytical prediction in Eq. ͑20͒ for the fixed flux boundary condition, that the system approaches a lower value for ͗B͘ ͑0.3 vs 0.4 for d = 24 and 0.7 for d =12͒, and correspondingly smaller value of d͗A͘ / dt.Also, ͗B͘ arrives in the neighborhood of this value much quicker than for shallower systems, but appears noisy and never asymptotes to a steady state.The deeper layer is also associated with larger, rapidly varying velocities, exhibiting a competition of unstable modes with new plumes appearing between already established plumes throughout the simulation.However, the plumes, both new and old, rapidly merge or migrate to the boundaries.In the latter stages of the simulation a large horizontal gradient in MB + emerges ͑in this simulation, mostly blue on the right and clear on the left͒, generating a single strong circulation cell which continuously shears several emerging plumes to such an extent that they are translated to the right and mixed by stretching and folding ͑see Fig. 10͒.
The initial dominant wavenumber is larger for shallow systems than in deeper ones, but in both cases the system ends up with one or two cells after the sequence of merging events.

A. Comparisons with experimental results, and linear and weakly nonlinear theory
Qualitatively, the results from simulations fit extremely well with experimental measurements of pattern initiation and wavenumber evolution ͓15͔.Common to both is an initial rapid increase of pattern wavenumber ͑usually a dou-bling͒ after the instability arises, coupled with initial thinning of the width of the blue downwelling plumes.Beyond this very early stage is a cascade of plume merging events which reduces the wavenumber and broadens the plumes.Overall, simulations and experiments reveal plume aggregation processes, which either draw the pattern to the center of the Petri dish and can lead to a localized pattern packet with slightly larger wavenumber, or lead to accumulations of plumes close to, or at, the boundaries.Furthermore, in deeper chambers large scale circulation is observed in simulations that shear the plumes.Such plume shearing is evident in many experiments, and can be observed quite clearly in Fig. 1͑b͒ as short radial line patterns when seen from above.
The simulations predict that for the standard parameter values and for a range of boundary conditions, the initial instability arises around 15-18 nondimensional units ͓see inset of Fig. 3͑d͒ for ͉͗v͉͔͘.This corresponds to a dimensional time of 3571-4286 s, which is in good agreement with the reference case experiment reported in ͓15͔ of 3100 s for pattern onset.Moreover, one can estimate the theoretical onset time from pseudo-steady state linear analysis in ͓20͔ as the time at which the characteristic growth rate approaches the natural time scale of the system, k obs .For the reference case in ͓15͔ this provides a pattern onset time of 3500 s.
Linear theory using the pseudo-steady state approximation ͓20͔ indicates that the dominant initial wavenumber should rapidly increase in the first moments after pattern onset.Simply put, this is due to an evolving neutral stability curve that possesses a zero most unstable mode ͑with zero growth rate at k =0͒.For a fixed value of the Rayleigh number, one may view the increasingly unstable GLA profile as inducing a neutral stability curve that everywhere drops and allows for a wider band of unstable modes with ever greater linear growth rates.This is observed both in experiments ͓15͔, where approximate doubling of dominant wavenumbers is common ͑although the three-dimensional experiments allow for more complicated increases͒, and in the twodimensional simulations, where wavenumber doubling is typically observed on a very fast time scale.Finally, we refer to the amplitude equation derived in ͓21͔ for a FFBC from a steady state, which took the form

͑23͒
where k 1 = −1.7,k 2 = −94, k 3 = 2.7, k 4 = 0.29, R c = 4.9ϫ 10 −4 ͑2 s.f.͒.Direct simulations of this equation ͓21͔ indicate that grain-coarsening ͑modes are unstable relative to modes of smaller wavenumber͒ is predominant after the onset of the instability, via a series of plume merging events, which take on a very similar appearance to the full simulation results for FFBC presented herein.As the amplitude equation has been derived from a steady state, which is generally far from being obtained in experiments or simulations, the initial increase in wavenumber is not observed.Neither is the horizontal interaction of plumes with the boundaries as the amplitude description is derived for containers of infinite horizontal extent.
The initial numerically predicted wavenumber is around k num = 0.02 nondimensional units, which equates to a physical wavenumber of 2k num = 0.13, or a wavelength of H ¯/ k num = 1.7 cm.This agrees with the experimentally observed initial wavelength of ͑12.4 cm/ 9.5=͒1.3cm ͓15͔.In considering this good agreement it also must be remembered that we are using independently measured parameter values with no fitting.The subsequent increase in wavenumber in simulation and experiment gives wavelengths of ͑0.031 cm/ 0.03 =͒1.03 cm and ͑12.4 cm/ 13.5=͒0.91cm, respectively.The long term simulation results reported in this paper are also consistent with the experimental measurements in ͓15͔.
Experiments indicate that oscillations are dominant in the long term in deeper containers ͓15͔.This is also our experience in the simulations.However, containers with aspect ratios such that d ϳ L or greater are beyond the scope of the current paper.

B. Minimization of "energy?"
It is tempting to make a rough "kinetic" and "potential energy" argument in line with the methods of ͓27͔ for a numerical study of bioconvection.Even though such descriptions are not rigorous they suggest a process of minimization based upon simulation measurements.Consider the definition of total "kinetic energy," K e , defined as half of the sum of local density times the square of velocity, such that where w n m represent Clenshaw-Curtis weights, and P e ͑"potential energy"͒ due to vertical gradients in GLA, or the "first moment" in the vertical direction of the deviation of gluconic acid with respect to its average, such that We measure the deviation of gluconic acid relative to its average ͗A͘, as the continuous growth of ͗A͘ hides the true useful potential energy, which is distributed inhomogeneously in space.
In Fig. 11 we plot K e vs P e for the standard parameter values with a range of oxygen boundary conditions.In all cases potential energy increases from zero to such a level that the instability is triggered and kinetic energy rapidly increases from zero to a high value, accompanied by a drop in P e .This is followed by a drop in K e and a more gradual increase in P e .Oscillations in these two energies follow, through a sequence of either circulation-circulation or circulation-boundary interactions that appear to generate attractors with smaller values of P e and larger K e .In general, this culminates in a spiral-type attracting fixed point for the final one-cell circulation pattern.This stable-focus-cascade structure is common to all of the simulation results.It suggests that given the number of circulation rolls, the system evolves towards a state with minimal total energy, but in the long term will evolve towards a simple reacting and convecting state that minimizes P e globally.

C. Impact on aggregate rates of reaction
It is clear that the impact on aggregate rates of reaction is governed by ͗B͘ as it is the rate of production of ͗A͘.In the simulations ͑e.g., Fig. 3͒ we see that initially the concentration profiles evolve such that ͗B͘Ϸ1, and ͗⍀͘ decreases from 1 until a point is reached where there is not enough oxygen to support ͗B͘Ϸ1.͗B͘ then drops resulting in lower rates of production of ͗A͘.The system then either approaches a horizontally homogeneous steady state, or becomes unstable before a steady state can be reached and develops a pattern.There are initial wiggles for all averaged quantities in the simulations exhibiting patterns.The upper oxygen boundary condition then determines how the aggregate rates of reaction are affected.
Beyond times with transient solutions for the FFBC, no impact of pattern generation is felt relative to the horizontally homogeneous steady state.This is due to the fact that a steady state in the average quantities ͗B͘ and ͗⍀͘ exists, as is revealed in Sec.V A and is the same as for the horizontally homogeneous steady state ͑i.e., ͗A͘ gets produced at the same rate͒.
In the regime before pattern onset, the FCBC case is similar to the FFBC case if the initial setup is such that the concentration and flux are matched to exhibit the same solution in the absence of pattern.However, the steady state analysis for averaged quantities does not follow as for FFBC and one can see wiggles in ͗B͘ as well as wiggles in the vertically averaged solutions in Fig. 7.
We see a drop in ͗B͘ ͑i.e., drop in production rate of ͗A͒͘ with the mixed boundary conditions for ␤ c = 1 and ␤ f =3 compared to the FFBC case and a greater drop for ␤ f =9.There is no significant change in production rate of ͗B͘ when ␤ f = 0.2 compared with the FFBC case.It should be noted again that ͗B͘ is more noisy for the mixed boundary conditions due to continual rebirth, shearing, and mixing of plumes.

VII. CONCLUSION
In this paper we have constructed a numerical procedure to simulate equations describing the instability mechanism for a chemically driven hydrodynamic instability.In particular, we have investigated pattern formation in the methyleneblue-glucose system and probed the underlying instability mechanism, namely the formation of dense regions of the fluid at the upper surface due to reactions generating the relatively but slightly dense product gluconic acid.We have considered the upper oxygen boundary condition in some detail in order to shed light on the discrepancy between previous linear stability analyses with fixed oxygen concentration upper boundary conditions ͑FCBCs͒ and weakly nonlinear analysis with fixed oxygen flux upper boundary conditions ͑FFBCs͒.Statistical measures constructed from the numerical solutions were used to aid discussion of the results, and displayed good agreement with theoretical predictions on averages from the equations, suggesting that the numerical procedure is accurate.Furthermore, we have performed a systematic numerical exploration of parameter space, and documented the extent to which solutions vary for realistic parameter values ͑noting that the reaction rates and diffusivities have been measured previously to a high degree of accuracy͒.
We find that solutions with a FFBC are very different to solutions with a FCBC.The former has plumes that proceed through a cascade of plume merging events.The plumes are typically horizontally symmetric but lack vertical symmetry, both for the plumes and in terms of downwelling and upwelling regions.In contrast, the latter generates plumes that typically do not merge and are repelled by one another ͑with occasional rapid crossing events of plume-antiplume dou-blets͒.These solutions also display much greater vertical symmetry ͑for individual plumes͒ but instead are horizontally asymmetric.Solutions with mixed boundary conditions share some of the features of each of the solutions.We note that the FFBC and mixed boundary conditions give qualitative solutions more akin to those seen in experiments, although it is clear that both are basic models of the real highly complex chemicophysical boundary condition.Imposing a FCBC appears to be a much harsher constraint, which adjusts flux so that oxygen depleted fluid never reaches the surface.
We find that the numerical solutions quantitatively agree with linear analysis ͓20͔ with a FCBC in terms of pattern wavelength, time scales, and critical Rayleigh numbers.Also, the predictions from a pseudo-steady linear analysis are consistent in that the wavelength of the pattern is seen to decrease initially for a short period of time after the instability.Furthermore, for the FFBC, the numerical solutions agree with weakly nonlinear analysis ͓21͔ in that for an intermediate period after the instability a grain coarsening cascade ensues through plume merging events.This case also agrees very well with what is observed experimentally ͓4,15͔ in terms of time and length scales of the pattern.However, it is clear that the linear predictions ͑FCBCs͒ and weakly nonlinear predictions ͑FFBCs͒ are at odds with one another.The answer is simply that there should be no agreement as the solutions with different boundary conditions are qualitatively distinct.
Of interest is that the form of the upper oxygen boundary condition could be measured indirectly by measuring any associated entrainment of oxygen in experiments.From the opposite perspective, controlling the precise form of the oxygen boundary condition could help to improve the yield of products in a chemical engineering setting.
A major extension of the work will be to consider threedimensional pattern formation in this system so that computed patterns can be compared with experiments.It is likely that such a numerical simulation would enable a more direct comparison with experiments and may suggest which form of boundary condition is more appropriate.This is work in progress.The methods in this paper should be extended to other chemical reactions and systems to determine the class of system for which solutions are highly boundary condition dependent, and to establish whether or not the form of solutions presented herein are generic.
FIG. 1. Experiments indicating a range of patterns, as seen from above.See ͓4,15͔ for more information and other examples.The first two patterns are in Petri dishes which are open at the upper surface ͑the images are 13.5 cm wide͒ and the third is in a completely filled container ͑of diameter 5.1 cm͒ with an oxygen permeable lid.
FIG.2.͑Color online͒ Five snapshots around the initial instability of simulations, for the standard parameter values with a fixed flux oxygen boundary condition at the upper surface ͑FFBC͒.Fluid velocity is superimposed upon each of the plots ͑if max͉u͉ Ͼ 10 −3 ͒, which are arranged in groupings of ͑a͒ A − ͗A͘, ͑b͒ B, and ͑c͒ ⍀.Snapshots are provided at times of 17, 18, 19, 20, and 21 nondimensional units from bottom to top.The maximum modulus of velocity is indicated on the right.Recall that there are 10 3 numerical time steps between each recorded integer time t d ͑Ϸ4 min͒.͑d͒ Plots of Var x ͑A͒ and Var x ͑B͒ ͑including data at time 1000 units͒ illustrating variations from a horizontally homogeneous solution ͑can be compared with perturbation solutions in linear and weakly nonlinear analyses͒.
FIG. 4. ͑Color online͒ Vertically averaged quantities ͑a͒ ͗A͘ z − ͗A͘, ͑b͒ ͗B͘ z , ͑c͒ ͗⍀͘ z , and ͑d͒ ͉͗v͉͘ z for the standard parameter values with a fixed flux oxygen boundary condition at the upper surface ͑FFBC͒.The plots on the left are for 0-1000 nondimensional time units, and on the right for 0-100 nondimensional time units.
FIG.6.͑Color online͒ Plume interaction and annihilation for the standard parameter values with a fixed flux oxygen boundary condition at the upper surface ͑FFBC͒.Snapshots are presented at nondimensional times of 120, 130 and 140.The images portray a plume merging event.
FIG. 7. ͑Color online͒ Results for standard parameter values with a fixed oxygen concentration upper boundary condition ͑FCBC͒.͑a͒ and ͑b͒ are the vertically averaged quantities ͗A͘ z − ͗A͘, ͗B͘ z , for early and long times, and ͑c͒ are snapshots of A − ͗A͘ ͑top͒ and B ͑bottom͒ at a time of 120 nondimensional units, with velocity superimposed.These images clearly indicate the horizontal asymmetry of cells.
FIG. 10. ͑Color online͒ Plumes developing in a strong shear flow arising in a deep layer of fluid, d =35 ͑70 collocation points in the z direction͒, for the standard parameter values with mixed boundary conditions ͑␤ c =1, ␤ f =3͒.The snapshot of ͗B͘ z is presented at a nondimensional time of 1000.

TABLE I .
Description of dimensional quantities, values, and units.
H ¯sublayer depth= ͱ D / k obs 0.031 cm direction ͑except where indicated otherwise͒, with a time step of 10 −3 t d .

TABLE II .
Nondimensional parameters, their definitions, and the standard parameter values to 3 s.f.

TABLE III .
Summary of general solution dependence on reaction and diffusion parameters for standard parameter values with mixed boundary conditions ͑␤ c =1, ␤ f =3͒.