In a gravitational field, a horizontal interface between two miscible fluids can be buoyantly unstable because of double diffusive effects or because of a Rayleigh-Taylor instability arising when a denser fluid lies on top of a less dense one. We show here both experimentally and theoretically that, besides such classical buoyancy-driven instabilities, a new mixed mode dynamics exists when these two instabilities act cooperatively. This is the case when the upper denser solution contains a solute A, which diffuses sufficiently faster than a solute B initially in the lower layer to yield non-monotonic density profiles after contact of the two solutions. We derive analytically the conditions for existence of this mixed mode in the (*R*, δ) parameter plane, where *R* is the buoyancy ratio between the two solutions and δ is the ratio of diffusion coefficient of the solutes. We find an excellent agreement of these theoretical predictions with experiments performed in Hele-Shaw cells and with numerical simulations.

## I. INTRODUCTION

Studies on the mixing of fluids due to buoyancy-driven convective instabilities are of generic interest because of the tantamount importance of such convective processes in technological and environmental problems. As an example of such instabilities, the well-known Rayleigh-Taylor (RT) instability occurs when a denser solution lies on top of a less dense one in the gravity field leading to interpenetrating fingers and vivid convective motions (Fig. 1(a)).^{1,2} RT convection is observed in applications as diverse as CO_{2} sequestration,^{3,4} stellar, and planetary interior dynamics,^{5} granular systems,^{6} or sediment and contamination transport in soils^{7} as well as in the atmosphere or oceans for instance. Other much studied buoyancy-driven instabilities are double diffusive instabilities observed in the presence of two adverse density gradients due to two different variables diffusing at different rates.^{8–10} Such differential diffusion phenomena are for instance at the origin of thermohaline convective motions in the oceans due to differential diffusivities of salt and heat^{11} and are also observed in magma crystallization phenomena^{12} or in stellar dynamics for instance. They can also occur as soon as two solutes with different diffusion coefficients (such as salt and sugar, for instance, see Figs. 1(b) and 1(c) participate in density gradients starting from an initially stable density stratification.^{13–17} Because of the large range of spatial scales involved in applications of RT and double diffusive instabilities, theoretical modeling backed up by laboratory-scale experiments have been developed to analyze their dynamics. Up to now, these classical buoyancy-driven instabilities have been characterized individually, but the dynamics resulting from the interaction between them has not yet been considered, in contrast to numerous other pattern-forming systems.

In this context, we propose a unifying classification of all such buoyancy-driven instabilities of a two-layer miscible system. We furthermore demonstrate both theoretically and experimentally the existence of a hitherto unexplored mixed mode (MM) of instability resulting from the coupling between RT and differential diffusion convective modes (Figs. 1(d)–1(f)). This mixed mode regime is observed in a stratification of a denser solution of a solute A on top of a less dense solution of another solute B provided A diffuses sufficiently faster than B. It is characterized by both a saturation of the deformation of the initial contact line in time and “Y” shaped antennas growing on both sides of it. We derive analytically the existence zone of this MM in the parameter space spanned by a buoyancy ratio *R* and a diffusion coefficient ratio δ. Excellent agreement between these predictions, experimental data, and numerical simulations of the relevant model are obtained.

## II. BUOYANCY-DRIVEN INSTABILITIES

We consider a two-dimensional vertical stratification of a solution of A at a concentration *A*_{0} and density ρ_{A} overlying a miscible solution of B at concentration *B*_{0} and density ρ_{B} with the gravity field pointing downwards along the *x* direction. The initial condition for the dimensionless concentrations *a* = *A*/*A*_{0}, *b* = *B*/*B*_{0} is taken as *a* = *H*(−*x*) and *b* = *H*(*x*), where *H* is the Heaviside function and *x* = 0 is the initial position of the horizontal contact line. The dimensional density ρ_{d} of the solution is assumed to vary linearly with the concentrations as ρ_{d} = ρ_{0}(1 + α_{A}*A* + α_{B}*B*) where the solutal expansion coefficients are

_{0}is the density of the pure solvent.

Classically, the convective instabilities affecting such a stratification of a fluid A on top of a miscible fluid B in the gravity field are classified in three regimes.^{10} The RT instability occurs when the density ρ_{A} = ρ_{0}(1 + α_{A}*A*_{0}) of the upper solution of A is larger than the density ρ_{B} = ρ_{0}(1 + α_{B}*B*_{0}) of the lower solution of B. The interface deforms in “fingers” of the less dense fluid rising between the denser sinking ones hence also the name “density fingering” for this instability, which has been much studied both theoretically and experimentally.^{1,2} Fig. 1(a) shows an example of a RT dynamics when a denser solution of glycerol is put on top of a less dense solution of sucrose in a Hele-Shaw cell, i.e., two glass plates separated here by a gap of 0.25 mm and oriented vertically in the gravity field.^{18}

If ρ_{A} < ρ_{B}, i.e., if we start from an initially statically stable stratification of a less dense fluid on top of a denser one, instabilities can nevertheless occur because of differential diffusion effects. These can take place if *D*_{A} and *D*_{B}, the diffusion coefficients of A and B are sufficiently different, i.e., if their ratio δ = *D*_{B}/*D*_{A} ≠ 1.

If δ > 1, the so-called double-diffusive (DD) instability (also commonly named “salt fingering” because of its general applicability in ocean dynamics^{8}) can be observed. The DD modes also produce ascending and descending fingers across the initial contact line^{10,13,14} as seen on Fig. 1(b), where a less dense solution of slow diffusing sucrose overlies a denser solution of faster diffusing glycerol.

Eventually, if ρ_{A} < ρ_{B} but δ < 1, diffusive layer convection (DLC) modes develop convective vortices at symmetric distances above and below the unperturbed interface as seen on Fig. 1(c) for a solution of the faster diffusing KCl above a denser solution of sucrose.^{10,16} DLC is due to the fact that *D*_{A} > *D*_{B} and hence, the component A in the less dense upper solution diffuses faster downwards than the component B of the denser lower solution diffuses upwards. This leads to both a depletion and accumulation zone developing, respectively, above and below the interface, and in which locally adverse density gradients trigger convection.^{10,15} The initial contact line remains unperturbed as seen on Fig. 1(c) as the DLC pattern develops symmetrically around it.

To unify the description of these various modes of convection, we non-dimensionalize the problem and classify the various convective regimes in a suitable parameter plane. In the absence of convection, the two solutions of A and B mix by diffusion and the base state of the problem is

*x*and

*t*the dimensionless space and time.

^{10}The corresponding dimensionless base state density profile is

_{d}− ρ

_{0})/ρ

_{0}α

_{A}

*A*

_{0}. The buoyancy ratio

*R*= α

_{B}

*B*

_{0}/(α

_{A}

*A*

_{0}), which is also equal to

*R*= (ρ

_{B}− ρ

_{0})/(ρ

_{A}− ρ

_{0}), expresses the relative contribution of species B and A to the dimensional density profile. In dimensionless units, the density ρ

_{a}of the upper solution is

_{b}=

*R*.

In the parameter space (*R*, δ), the RT instability occurs when ρ_{a} > ρ_{b}, i.e., in the half plane *R* < 1 (Fig. 2).^{10} The DD and DLC modes can exist in the reverse situation *R* > 1 in the δ > 1 and δ < 1 regions, respectively. Stability regions and dynamical properties of these three RT, DD, and DLC regimes have already been described previously^{10} and we focus here only on the lower half plane δ < 1.

To further investigate the influence of differential diffusive effects on the dynamics when δ ≠ 1, we note that the density profile features non-monotonic behavior^{19–21} for δ > max(1, *R*^{2}) (not shown here but detailed in Ref. 10) or δ < min(1, *R*^{2}). In particular, if δ < *R*^{2} < 1 (Fig. 2), a maximum ρ_{max} and minimum ρ_{min} in density develop locally at η = ±η_{m} where

As seen on the sketches of Figs. 2 and 3, at fixed δ < *R*^{2} < 1, increasing *R* simultaneously decreases the density difference Δρ_{0} = ρ_{a} − ρ_{b} = 1 − *R* between the two pure solutions and increases Δρ′ = ρ_{max} − ρ_{min}, the density difference between the two extrema. Let *R*_{c} denote the critical value of *R* when Δρ_{0} = Δρ′. Fig. 3(a) shows that the local maximum in the base state density profile has then a peak equal to 1 and a local minimum equal to *R*, the blue dashed line representing the initial density profile. The location η_{m} at which ρ_{max}(η_{m}) = 1 or also

*R*

_{c}→ 0.5 when δ → 0 while

*R*

_{c}→ 1 as δ → 1. An expansion about δ = 1 by writing

*R*

_{c}= 1 − ε and δ = 1 − ε

*n*and letting ε tend to zero, yields, to first order in ε, an equation for

*n*. By writing

*n*= 2/(1 − 2

*m*

^{2}), we obtain to leading order the equation

*n*≈ 6.8. One notes that the equation δ =

*R*

^{n}is consistent with this expansion near

*R*= δ = 1. Fig. 2 delimits the regions in the (

*R*, δ) parameter plane where the various types of density profiles are obtained for δ < 1.

## III. MIXED MODE

To link the above theoretical density profiles to experimental observations and test the importance of the non-monotonic character of the density profile on the nonlinear dynamics, we have performed experiments in Hele-Shaw cells scanning various values of δ and *R* (Figs. 1 and 2). To analyze the dynamics at different values of δ, we use KCl above sucrose (δ = 0.27), NaCl on top of sucrose (δ = 0.34), KCl above glycerol (δ = 0.53), and glycerol overlying sucrose (δ = 0.51). The diffusion coefficient and solutal expansion coefficient of the solutes used in the experiments are summarized in Table I. To vary *R* for one specific pair of solutes (fixed δ), we use initial concentrations *A*_{0} and *B*_{0} in the order of 10^{−2} mol/l. The dimensional density ρ_{A,B} of these solutions is measured using an Anton-Paar densimeter. One of the two solutions is then progressively diluted and its density measured to vary *R* according to the formula *R* = (ρ_{B} − ρ_{0})/(ρ_{A} − ρ_{0}) where ρ_{0} is the density of water. This procedure ensures reproducible experiments independently of slight changes of density with temperature and pressure.

As color indicator or dyes can perturb the dynamics,^{25,26} visualization is made by a Schlieren technique^{27} tracking changes in the index of refraction induced by the dynamics. In effect, the experimental images (Fig. 1) show changes in the gradient of the refraction index in a grey scale ranging from white to black between the absolute minimum and maximum values reached during the dynamics respectively. Initially, the contact line corresponds to the maximum black value in this scale as it is the location of the largest gradient. We follow the dynamics of this “contact line” at later time by following the evolution of this maximum.

In the DLC regime, i.e. if *R* > 1 for δ < 1 (black squares in Fig. 2), the local minimum above and local maximum below the contact line for an initially stable density stratification (Δρ_{0} < 0) lead to the dynamics of Fig. 1(c). The contact line is initially non-deformed and convection develops at symmetric distances above and below it in the zones of the locally statically unstable density gradients.

In the RT zone for *R* < 1, the dynamics depends on the value of δ. For δ > *R*^{2}, the density profile ρ is monotonically decreasing along gravity (Fig. 2), and a classical RT dynamics is observed, with regular smooth fingers and a sinusoïdal deformation of the initial contact line at the onset of the instability such as in the case of a denser solution of A on top of the pure solvent studied experimentally by Fernandez *et al.* in a Hele-Shaw for instance.^{2}

If δ < *R*^{2}, the density profile becomes non-monotonic (Fig. 2) and differential diffusion effects start to impact the RT instability. As long as *R*_{c} < δ < *R*^{2}, then Δρ_{0} > Δρ′, i.e., the unfavorable density difference between the two pure solutions of A and B is larger then the locally stable zone in the middle. Hence, the RT mode dominates and strongly deforms the contact line. Differential diffusion effects manifest themselves only by deforming the upward moving caps of the fingers in a mushroom way (see Fig. 1(a)). This is reminiscent of what is observed in nonlinear simulations of RT instabilities with a fast diffusing upper component performed by Trevelyan *et al.*^{10} even though the refractive index images of Fig. 1 cannot directly be compared to their numerical concentration maps. In this regime, the mixing length defined as the distance between the furthest upward and downward location of the interface deformation grows continuously in time (Fig. 4).

If δ < *R*_{c}, then Δρ_{0} < Δρ′ such that the values of the extrema of the non-monotonic density profile are outside the range of the end-point values of densities for the pure fluids (see Fig. 3(b)). The dynamics features then a combination of RT and DLC trends giving rise to the mixed mode. The DLC characteristic leads to little “Y shaped” convective structures developing around the interface at the locations of the local adverse density gradients between the end-point values and the extrema (Figs. 1(d) and 1(e)). These “antennas,” which keep on growing in time resulting in a steadily increasing mixing length (Fig. 4(b)), alternate on average upwards and downwards out of the interface, which initially remains almost unperturbed. Later on, the RT mode comes into play as well due to the fact that ρ_{a} > ρ_{b} and the interface itself starts to be modulated (Fig. 1(f)). However, the amplitude of this modulation saturates after a given time as seen on Fig. 4(c). This saturation is different from the continuous growth characteristic of the dominance of the RT mode for the cases where δ > *R*_{c}.^{10} The red diamonds in Fig. 2 represent the experimentally determined MM regimes where the DLC contribution to the dynamics begins to dominate over the RT instability leading to a saturation of the amplitude of the interface deformation after a given transient. At a fixed δ < 1, the contribution of the RT mode on the dynamics decreases as *R* increases and δ decreases. The related smooth decrease in the mixing length from RT to the DLC mode agrees with the numerical trends computed numerically in Trevelyan *et al.*^{10} (see Figure 18(b) in Ref. 10). This is also the case for the saturated amplitude of the interface modulation, which decreases when *R* increases as seen on Fig. 4(c).

The wavelength of the mixed mode, measured as the distance on a horizontal line between two consecutive upside (or downside) “Y antennas” does not remain constant in time mainly due to the waving of the interface, which induces a motion and sometimes even a possible collision of the plumes. The wavelength depends on the (*R*, δ) values and is roughly equal to 1.3 mm for the experiment shown in Figs. 1(d)–1(f). For a given couple of solutes, i.e., for a given value of δ, the experimental wavelength slightly increases upon increasing R.

In short, differential diffusion effects can deform RT dynamics when *R* < 1. As long as the extrema of the non-monotonic density profiles are not larger than the end-point densities of the pure solutions (*R*_{c} < δ < *R*^{2}, green circles of Fig. 2), the RT mode remains dominant, the interface strongly deforms and the mixing length nearly equals the amplitude of the interface modulation and both keep growing in time. In the MM regime (δ < *R*_{c}, red diamonds of Fig. 2), these extrema are outside the end-point values and the local intermediate stabilizing zone stops the continuous growth of the interface. This induces a saturation of the growth of the interface modulation with DLC antennas growing on both sides where the local adverse gradients operate, leading to a mixing length still increasing but at a rate decreasing when δ decreases or *R* increases.

## IV. NONLINEAR SIMULATIONS

Signatures of the non-monotonic density profile and of the MM mode can been obtained in nonlinear simulations of a diffusion-convection model of the problem as well. In the modeling, we assume that the solutions are dilute so that the Boussinesq approximation can be made and the fluid can be treated as incompressible. The evolution equations are the continuity equation, an evolution equation for the flow velocity

*a*,

*b*of the solutes A and B. To simulate the above experiments performed in Hele-Shaw cells, we use two-dimensional Darcy's law for the flow equation but a generalization to the Navier-Stokes equation could be done along the same line. In the absence of a given physical length scale, the problem is non-dimensionalized using the characteristic velocity

*u*

_{c}=

*gK*α

_{A}

*A*

_{o}/ν, where

*g*is the magnitude of the gravitational acceleration,

*K*the permeability, and ν the kinematic viscosity defined as μ/ρ

_{0}with μ being the dynamic viscosity. Then the characteristic length,

*l*

_{c}, and time,

*t*

_{c}, are given by

*D*

_{A}/

*u*

_{c}and

^{10}

^{,}

^{28}In the numerics, the DLC mode has a clear signature: the stream function features two rows of alternating vortices around the initial contact line instead of one single roll structure across the interface characteristics of the RT instability.

^{10}The transition from a dominant RT instability to MM is shown on Figure 5 with the concentration map of species B for R = 0.85 and δ = 0.8, 0.6, 0.4, and 0.2 from left to right, respectively. The colormap has been chosen in order to track the isoconcentration 0.5 as a marker of the evolution of the initial interface. For δ = 0.8, i.e., δ >

*R*

^{2}, the density profile is monotonic decreasing along the gravity field and a RT mode strongly deforms the interface right from the beginning. For δ = 0.6 and 0.4 (green circles on Fig. 2), the coupling between RT and DLC modes provides rising plumes and sinking fingers developing around a saturated interface deformation, the amplitude of which decreases with δ. The last panel for δ = 0.2 (in the MM red diamonds zone below

*R*

_{c}on Fig. 2) shows the mixed mode dynamics with a small RT deformation of the interface of asymptotic constant amplitude and DLC convective motions developing around it where the two extrema of the non-monotonic density profile are located. The evolution of the mixing lengths (computed in Ref. 10) and of the deformation of the interface for a constant δ and increasing

*R*is similar to those obtained in experiments. In Fig. 5, the widths and times of images have been chosen to highlight the scales on which the patterns have similar longitudinal extensions. This shows that, when the RT modes largely dominate the dynamics (δ = 0.8 and 0.6), the patterns grow on similar time and space scales. On the contrary, when the DLC mode starts to dominate (δ = 0.4 and 0.2), the wavelength of the pattern decreases and a similar extension is reached quicker, in other words the system is more unstable when δ decreases.

Note that the existence of the “Y” antennas is more easily obtained in reconstructed density profiles (Fig. 6) than in concentrations maps. To understand the origin of these “Y” shaped protrusions, consider an interface perturbed by some sinusoïdal perturbation. The faster diffusion of denser A downwards is not compensated fast enough by the slower diffusion of the less dense B upwards, which creates denser minima below the interface and less dense maxima above it (Fig. 6 (top left)). Within the minima, A keeps diffusing out perpendicularly to the surface of the protrusion quicker than B can fill it, which progressively dilutes the inner part of the minimum, accumulating A on the outer rim of the interface (Fig. 6 (bottom left)). However, at the tip of the minimum, because of curvature effects, the inflow of B is larger and the outflow of A is smaller than at the sides such that an averaged density is obtained there. As a consequence, the denser lateral zones sink on the sides of the less dense rising middle zone, which gives the characteristic “Y” shape (Fig. 6 (right)). A symmetrical argument can be applied to the upper growing maxima.

## V. CONCLUSIONS

We have provided experimental evidences of the influence of differential diffusion effects on the RT instability when a denser solution of a fast diffusing species A is put above a less dense solution of a slower diffusing component B. If *R*_{c} < δ < *R*^{2} < 1, the RT fingers feature mushroom caps while if δ < *R*_{c} < 1, a mixed-mode instability combining properties of the RT and DLC modes is observed. The MM dynamics features “Y shaped” convective structures characteristic of the DLC mode growing both above and below the initial contact line modulated with a finite saturated amplitude. In the nonlinear dynamics, the mixing length giving the total extent of the fingered zone increases in both RT, MM, and DLC regimes with a slope decreasing as *R* increases. The amplitude of the contact line deformation on the other hand increases in the RT mode, saturates to a constant value in the MM regime while is vanishing in the DLC mode. These characteristics are recovered in numerical simulations as well.

The discovery and experimental characterization in 2D Hele-Shaw cells of the new mixed mode buoyancy-driven instability calls for future studies including simulations of the MM dynamics using 2D and 3D Navier-Stokes equation and experimental studies of it in 3D tanks. Search for analogous mixed modes in dynamics of reactive systems^{15,29} or in other hydrodynamic instabilities^{30} where differential diffusion effects have been evidenced would also be of particular interest.

Eventually our classification of all possible buoyancy-driven instabilities of a two-layer system in a simple (*R*, δ) parameter plane paves the way for a unification of future comparisons between experiments and simulations even in systems with very different length and time scales.

## ACKNOWLEDGMENTS

A.D. acknowledges Prodex, the “Actions de Recherches Concertées CONVINCE” programme and FRS-FNRS for financial support. J.C.-L. thanks MICINN for funding through research (Project No. FIS2010-21023) and the FPI grant associated to FIS2007-64698.