We performed an experimental and numerical investigation of a convective buoyancy-driven instability that arises during the injection of a denser miscible fluid into a less dense one in a rectilinear geometry. We visualized the instability using a shadowgraph technique, and we obtained quantitative information using micro-Particle Image Velocimetry. Numerical simulations provided further insights into the three-dimensional (3D) velocity field. We have shown that the instability only occurs above a certain Péclet number, Pe, depending on the Rayleigh, Ra, and Schmidt, Sc, numbers. We suggest scalings of the critical time, TC, and dimensionless wavelength, λ / h, of the instability, both of which increase with increasing Pe and Ra. Finally, we investigated the interactions of the instability vortices with each other and the geometry boundaries.

Hydrodynamic instabilities are a classic, yet constantly active and fascinating area of research in fluid dynamics. They appear in a wide variety of applications such as porous media mass transfer,1 geochemical flows,2 CO2 sequestration,3 H2 storage,4 or Li-ion battery design.5 Hydrodynamic instabilities can also be part of a chaotic mixing mechanism in microfluidics, that is, in the form of viscous fingering or alternating injection of fluids,6,7 and play a critical role in certain micro-separation processes.8,9 Recently, a trend of investigating buoyancy-driven phenomena on Earth and their absence in space is also dynamically developing, in view of the future possibilities of transferring relevant processes in orbit or in outer space (i.e., during long-term space travel).10,11 In this context, new techniques for material synthesis in the microscale emerge in the absence of buoyancy phenomena.12 

It has been observed that when a fluid 1 of density ρ1 and viscosity μ1 is injected into a miscible fluid 2 of density ρ2 and viscosity μ2, a tongue-like profile is created in the contact zone of the two fluids, where an abundance of instabilities may occur13–15 whether they are convective,16–18 viscous,19,20 wetting-front,21 or even Marangoni-driven.22,23 When this displacement takes place in a gravity field and μ 1 μ 2, then, a buoyancy-driven convective instability dominates the dynamics of the system.16,24 Generally, when a denser fluid is superimposed on top of a lighter fluid, the horizontal interface is gravitationally unstable, resulting in the well-known Rayleigh–Taylor instability.25,26 In the presence of a displacement flow perpendicular to the gravity vector, advection and buoyancy synergize toward the generation of a more complex instability27 as depicted in Fig. 1.

FIG. 1.

Schematic representation of the miscible displacement and the subsequent emerging instability mechanism: lateral view (a) and transversal view (b).

FIG. 1.

Schematic representation of the miscible displacement and the subsequent emerging instability mechanism: lateral view (a) and transversal view (b).

Close modal

Garik et al.28 first attempted to interpret this phenomenon as a Marangoni-induced instability, i.e., driven by interfacial tension gradients, while studying cellular mixing. Lister and Kerr29 and Didden and Maxworthy30 in their initial experiments attributed the instability forming at the tip of the gravity current to dust contamination. A splitting instability caused by buoyancy effects was observed in capillary tubes,31 but due to spatial confinement, the phenomenon did not produce more than one set of instability rolls in the flow cell. Obernauer32 first conducted experiments in radial and linear larger-scale geometries and suggested that the instability depends on the density difference between the two miscible solutions, the flow rate and the tilting angle of the flow cell with respect to the perpendicular, confirming the gravitational nature of the instability. A similar effect was observed during the formation of precipitate particles33 in the contact zone of two reactant solutions during radial displacement. Patterns in confined flow-driven microemulsion generation34 have been likewise attributed to such convective instabilities. Later on, Haudin et al.16 in their work concluded that the phenomenon is, indeed, buoyancy-driven and that the viscosity difference between the two fluids that participate does not affect the emergence of the instability (for the viscously stable displacement case, that is when the more viscous fluid is injected into the less viscous one). This conclusion was further tested and confirmed in microgravity experiments.35 In continuation of this work,16 Pótári et al.36 described the phenomenon analytically and coupled it to experimental observations present in precipitation patterns. The same system has been revisited by Balog et al.37 Subsequently, scaling laws for high Pe number cases have been proposed.36 More recent research work investigated the stabilizing effect of the walls in such geometries without any displacing flow. It has been pointed out that there is a possibility to stabilize an otherwise unstable stratification by sufficiently decreasing the vertical lengthscale of the system.38 Numerical simulations by Talon et al.39 and John et al.40 showed the existence of a buoyancy-driven instability in single viscous fingers. However, the effect of viscosity gradients on the generation of the instability was not completely separated from that of the buoyancy contribution in this system.

On this basis, to gain further knowledge on how this instability evolves in a narrow rectilinear cell, we studied the viscously stable displacement of a denser fluid into a less dense and fully miscible one in rectilinear geometries. We conducted experiments in a horizontally oriented flow cell with varying gapwidths, flow rates, and differences in density between the two fluids. Detailed measurements with a shadowgraph optics (Sec. III A) allow us to eliminate any possible side effect of dyes, as used in the previous studies16,32 for flow visualization. We characterize in detail the special type of patterns that emerge in the said unstable displacements (Secs. III A and IV C), and we introduce the corresponding spatiotemporal scalings (Sec. IV B). Additional numerical simulations on selected points in the experimental parameter space, based on the experimental observations, support our findings. The velocity fields from micro-Particle Image Velocimetry and numerical simulations provide a detailed understanding of the flow structure and the temporal evolution of the instability (Sec. IV D).

As depicted in Fig. 1, we consider a small flow geometry of height, h, and width, w. The geometry is initially filled with fluid 1 with density ρ1. For t = 0, the injection of fluid 2, with density ρ 2 > ρ 1, starts with a constant flow rate, q. In the absence of buoyancy effects, the injection flow results in the creation of a Poiseuille-like velocity profile, with average flow velocity, U = q / h w. The subsequent stretching of the border between the two miscible fluids gives rise to a gravitationally unstable stratification in the horizontally placed flow cell, as shown in Fig. 1(b).

As previously described,16 the unstable stratification in combination with the advection on the y-direction causes the appearance of multiple streamwise counter-rotating pairs of convective rolls. The general mechanism of the creation of the longitudinal streamwise-oriented rolls is well described in the literature.16,39,41 The counter-rotating rolls can be directly compared to similar convective vortices created for the sheared Rayleigh–Benard convection case.42 

To formulate the problem mathematically, we consider the case where a more dense fluid displaces a less dense fluid in a rectilinear geometry, as shown in Fig. 1. The two fluids are miscible in any proportion, and their viscosities are considered equal. Based on the approach of various similar studies of the instability,43,44 the phenomena are governed by the following three-dimensional, incompressible Navier–Stokes equations:
· u = 0 ,
(1)
ρ ( u t + u u ) = μ 2 u p + ρ g ,
(2)
where ρ is the density, μ is the viscosity, and g is the gravitational acceleration.
In combination with the convection-diffusion equation for the species that controls density, we obtain
c t + u · c = D 2 c
(3)
with
ρ = ρ 1 + c Δ ρ ,
(4)
where c = 0 1, and
Δ ρ = ρ 2 ρ 1 ,
(5)
where D is the diffusion coefficient of the density-controlling species, while ρ1 and ρ2 are the density values of the pure species 1 and 2, respectively.
Considering the hydrostatic and dynamic pressure contribution, we get
p = p + ρ g z ,
(6)
where g = g e z ̂.

Let us then define the characteristic scales as L * = h , T * = h 2 / D , U * = L * / T *, and P * = D μ / h 2, and the dimensionless quantities ζ, t ̃ , u ̃ , p ̃.

We then get
· u ̃ = 0.
(7)
By taking advantage of the Boussinesq approximation, we get from Eq. (2)
1 S c ( u ̃ t ̃ + u ̃ u ̃ ) = 2 u ̃ p ̃ R a ζ c ζ
(8)
and from Eq. (3)
c t ̃ + u ̃ · c = 2 c
(9)
obtaining the dimensionless parameters in the form of the Schmidt number, S c = μ / ρ D, and the Rayleigh number, R a = Δ ρ g h 3 / μ D. Furthermore, the nondimensional velocity takes the form of the Péclet number, P e = h u / D. For the experiments, the density and viscosity values of the denser solution in each case were used to determine Sc and Ra.

We have conducted experiments in a horizontally placed rectilinear geometry shown in Fig. 2. In contrast to previous experimental or numerical studies,24,32,36 where the geometry was considered nearly 2D, in our case, we have used Hele–Shaw cells with smaller h/w ratios. This fact allows for the development of a parabola-like velocity profile in the y-direction. Therefore, we also consider this aspect, which is present to some extent in previous experiments26,32 but assumed to be negligible. The flow cell comprised two optical glass plates separated by a thin polytetrafluoroethylene (PTFE) sheet of variable thickness, h = 0.12, 0.25, 0.5, 0.8, 0.9, 1.0  mm (more information on the dimensions of the flow cells are included in the SM). This setup allows for vertical optical access. The cell is first completely filled with one fluid; afterward, the second fluid is injected. The detailed values of the experimental and numerical trials are included in the SM. The cell is equipped with inlet and outlet ports for fluid handling (not shown in Fig. 2 for simplicity). The fluid is injected through an injection valve that is connected to a syringe pump (PHD ULTRA™, Harvard Apparatus, Holliston, MA, USA) with the use of PTFE tubing material. A wide range of flow rates (q) was offered by the syringe pump. To guarantee full miscibility and control the density of the injected solution, we used aqueous solutions with variable concentration of glycerol as solute. The used solutions along with their corresponding physical properties at 25 °C are included in Table I.

FIG. 2.

Sketch of the experimental setup used in this study: shadowgraph setup (a) and detail of the experimental Hele–Shaw cell (b).

FIG. 2.

Sketch of the experimental setup used in this study: shadowgraph setup (a) and detail of the experimental Hele–Shaw cell (b).

Close modal
TABLE I.

Physical properties of fluids used.

Fluid Glycerol (% v/v) Density ( g / cm 3) Viscosity ( cP)
0.9972  0.92 
6.8  1.0182  1.2209 
13.6  1.0363  1.5445 
Fluid Glycerol (% v/v) Density ( g / cm 3) Viscosity ( cP)
0.9972  0.92 
6.8  1.0182  1.2209 
13.6  1.0363  1.5445 

The density and viscosity values were obtained using an Anton Paar SVM 3001 pycnometer-viscosimeter (Anton Paar GmbH, Graz, Austria). To avoid any influence of a potential viscous-fingering instability, the flow cell was always pre-filled with the less viscous liquid (i.e., distilled water) and the more viscous fluid (glycerol solution) was injected. As the difference in the viscosity values among the used solutions is moderately small, the viscosity and density values of the denser solution were used to calculate all non-dimensional quantities (i.e., Sc, Ra) for each displacement combination.

To visualize the instability, we employed a shadowgraphy method which is sensitive to the second derivative of the refractive index field (and thereby sensitive to the second derivative of the h-averaged concentration field).45 We made use of a custom-made shadowgraph device (TSI Optics, Pulsnitz, Germany) in combination with a JAI GO-5100M, 2464 × 2056 px CMOS camera (JAI A/S, Copenhagen, Denmark).

In addition, we performed preliminary velocity field measurements employing a micro-Particle Image Velocimetry (μ-PIV) setup. The study was conducted using the same setup and parameters as for the shadowgraph experiments. The integrated μ-PIV system (LaVision GmbH, Göttingen, DE) illuminates the flow geometry with the use of an ND:YAG laser (Photonics Industries International Inc., Ronkonkoma, NY, USA). Fluorescent particles, with a mean diameter of d p = 5.03 μ m (Microparticles GmbH, Berlin, DE), dispersed in water and glycerol solutions were excited by the laser's pulse (532  nm) and re-emitted red light with a wavelength of 607  nm. The light emitted by the particles passes through a filter where all scattered background signal is removed. After magnification by a confocal microscope with a 3.2× objective lens (Zeiss Stereo Discovery V8, Carl Zeiss Microscopy GmbH, Jena, DE), the image is captured by the CMOS camera with a resolution of 1280 × 800 px (Phantom VEO 410L, Vision Research, Inc., NJ, USA). The image capturing and acquisition was synchronized and controlled by the DaVis software (LaVision GmbH, Göttingen, DE). The applied optics setup resulted in a depth of field of 0.18 mm. Hence, the velocity gradients in the vertical (z) direction are expected to influence the PIV results so that the measurement data mainly provide the approximate order of magnitude of the velocity field. The image acquisition rate was adjusted depending on the flow rate q and ranged between 0.2 and 1 kHz.

Computational Fluid Dynamics simulations were performed for a rectangular domain corresponding to the experimental setup (cf. Fig. S1 in supplementary material), using the twoLiquidMixingFOAM solver of the OpenFoam software package, a Volume-Of-Fluid solver. This solver is designed for modeling two miscible incompressible fluids. The incompressible Navier–Stokes equation [Eq. (2)] with the Boussinesq approximation is solved using the PISO method coupled with the scalar transport equation [Eq. (3)]. From the calculated concentration field (c), the density (ρ) can be corrected by applying Eqs. (4) and (5). The same method is used for the calculation of viscosity. Additionally, for all computations, the diffusion coefficient, D, was set to 10 9 m 2 s 1 for the species that controls density, in our case, glycerol. For temporal integration, we used the Euler method with a time step in the range of Δ t = 0.5 2 ms depending on the initial velocity. The time step was adapted to account for the CFL criterion, with the Courant number kept below 0.5. At the inlet and outlet, a constant initial velocity was set, and on the walls, a no-slip boundary condition was enforced. With regard to the concentration field at t = 0, a concentration of c = 0 was maintained for the outlet and internal space, c = 1 for inlet, and no flux boundary condition for the walls. The rectangular geometry was divided into hexahedral finite-volume cells with a spatial discretization of Δ x = 0.2 mm , Δ y = 0.05 mm , and Δ z = 0.04 mm. A grid convergence study was performed with a refined mesh (50% higher resolution for Δ x , Δ y, and Δ z, and a time step, Δ t = 0.5 s), which showed matching computational results. In order to directly compare the experimental shadowgraph images with the numerical simulation results, a “synthetic” shadowgraph image is reproduced from the simulation results of the corresponding case, where the magnitude of the 2 c z field is shown, with cz being the z-averaged concentration values.

In Fig. 3, shadowgraph experiments of a 6.8 % glycerol solution injected into distilled water in a flow cell ( h = 0.8 mm) are presented for different Pe numbers (i.e., injection flow rates). Under the shadowgraph, zones with a steep change in the refractive index second spatial derivative appear significantly darker or brighter than the background, allowing us to visualize structures in the concentration field during the displacement. For the lowest Pe experiment [Fig. 3(a)], no visual structures are observed, as no instability is present. As clearly seen in Figs. 3(b) and 3(c), the instability sets in and the counterrotating vortices can be visualized as stripes with brighter and darker areas than the background, following the displacement from left to right. The visual structures and the interaction of the vortices formed due to the instability will be further discussed in Secs. IV C and IV D. A preceding parabolic-shaped optical disturbance because of the Poiseuille profile of the displacement is visible in Fig. 3(c) (pointed with red arrows). This feature appears clearly in that case because of the high density contrast between the two solutions which is conserved for longer times due to the higher Pe. The parabolic profile also affects the evolution of the instability vortices, i.e., the vortices toward the middle of the cell are advancing quicker than the ones near the walls. For Figs. 3(a) and 3(b), the relative progress of the displacement can be estimated by the brighter areas near the walls [red arrows in Fig. 3(a)]. The no-slip condition causes a long-lasting concentration gradient near the walls, which is only smeared out later due to diffusive mass transfer.

FIG. 3.

Shadowgraph images of three experiments with Ra = 82 220 and increasing Pe: (a) Pe = 144, (b) Pe = 202, and (c) Pe = 865. The images are all taken at  10  mm downstream the inlet. Red arrows indicate the gradient near the walls in (a) and the Poiseuille profile preceding the instability in (c).

FIG. 3.

Shadowgraph images of three experiments with Ra = 82 220 and increasing Pe: (a) Pe = 144, (b) Pe = 202, and (c) Pe = 865. The images are all taken at  10  mm downstream the inlet. Red arrows indicate the gradient near the walls in (a) and the Poiseuille profile preceding the instability in (c).

Close modal

In the following, the relationship between the emergence of the instability and the flow conditions (Pe) is outlined and compared with preexisting research. It is observed that the higher the Pe is, either with an increase in h and/or q in our experiments, the faster the instability appears which leads to the creation of more vortices for the same flow cell width w (the wavelength decreases). When observing experiments under the same flow conditions (same Pe) but with different Δ ρ between the two fluids, the influence of density gradient is apparent (Fig. 4). Larger Δ ρ leads to an earlier onset of the instability and an increase of the number of counter-rotating vortices in the flow domain. Hence, the emergence and development mechanism of the instability clearly depends on the magnitude of buoyant forces as well.

FIG. 4.

Temporal progression of the instability through shadowgraph images of two experiments with Pe = 342 with (a) Ra = 30 260 and (b) Ra = 20 070. The images are all taken at  10  mm dowstream the inlet.

FIG. 4.

Temporal progression of the instability through shadowgraph images of two experiments with Pe = 342 with (a) Ra = 30 260 and (b) Ra = 20 070. The images are all taken at  10  mm dowstream the inlet.

Close modal

A similar trend for the instability wavelength depending on h has been previously observed for experiments with radial injection16 and numerical simulations in a rectilinear geometry with free slip walls.36 Also, note the faster progression of the displacement on the left side of Fig. 4. Even though the average flow velocity, U, and the Pe number are equal for both cases, the higher Δ ρ affects the displacement progress. This probably can be ascribed to a stronger gravity current occurring in that case. As discussed by John et al.,40 such a gravity current leads to a faster tip advancement. This fact might also contribute to the faster onset of the buoyant instability because steeper local density gradients are created. An analysis based on an advancing gravity current has been discussed by Pótári et al.36 

Using the introduced non-dimensional parameters, we can construct a flow map for the linear displacement as shown in Fig. 5(a). The Pe vs ScRa diagram is based on the collective data through the complete range of Δ ρ, h, and q. The slight changes in Sc are caused only by the variation of the glycerol concentration in the solutions used. The ScRa quantity is the product of the Sc and Ra numbers, also expressed as ScRa = Δ ρ g h 3 / ρ D 2. In the chart, ° indicates the cases where no instability was observed during the displacement, while + is used for the cases where the instability is present. The results of the numerical simulations, which are in good agreement with the experimental results, are also plotted in Fig. 5(a). In Fig. 5(b), the experimental correlation between Pe and Ra for the critical values is sketched. The cases with the lowest Pe in which instability occurs under otherwise identical experimental conditions (h and Δ ρ) are coined critical. For high ScRa values, which means that the magnitude of buoyancy effects is stronger, lower Pe values are required to destabilize the displacement. After fitting these critical values, we observe a correlation between the two of the form: P e ( ScRa ) 0.39. We can, thus, assume that the critical condition in our experiments is mainly controlled by the characteristic length (h), the density difference ( Δ ρ), and the magnitude of the advective flow (U). All of the above can be combined using the quantities expressed by ScRa and Pe.

FIG. 5.

Flow regime chart for the displacement including experimental and numerical values (a) and correlation of the experimental critical points for the instability emergence (b).

FIG. 5.

Flow regime chart for the displacement including experimental and numerical values (a) and correlation of the experimental critical points for the instability emergence (b).

Close modal

For experiments at low h values (0.12 and 0.25 mm), we observed no instability at any flow rate. This is seemingly an effect of spatial suppression of the instability (small h), as it is the case for a Rayleigh–Taylor convection observed in similar miscible systems without advection.38 In that case, the diffusion timescale is smaller than the instability onset time, smearing out the glycerol concentration gradient, and as a result, the density gradient as well. Because of this confinement stabilization effect, the prediction of the critical Pe number, where the instability arises, should not be extrapolated for lower ScRa values, as such instability might not arise for ScRa values below a threshold limit.

The onset time, tc, of the instability, defined as the time interval between the entrance of the injected fluid in the flow cell and the moment the instability sets in, seems to vary when the experimental conditions change. To draw a more consistent comparison, a non-dimensional onset time is defined as T c = t c / T *. To define the physical onset time, tc, the standard deviation of the grayscale value of the shadowgraph images (or the 2 c z, for the CFD results) is used. The tc is approximately set to the time when the standard deviation starts increasing rapidly (cf. see supplementary material).

In Fig. 6(a), the instability onset time is plotted against the Pe number, for various experimental conditions. Each data series represents experiments carried out for the same h and Δ ρ, while the only varying parameter affecting Pe was the flow rate, q. Tc shows a weak proportional decrease according to 1 / P e, a finding in accordance with previous studies in similar systems.24 For systems without advection, scaling laws of similar form have been proposed: either for the stationary Rayleigh–Taylor instability38 or for a solutal system, based on the Rayleigh–Bénard scaling laws.46,47 Although other studies suggest an independency of tc by Pe,16,36 we suppose this is caused by the relatively high flow rates used in those cases, where the differences between tc might be negligible.

FIG. 6.

Dependence of the non-dimensional critical time, Tc on Pe for various experiments with different Ra (a). Correlation of Tc with Pe and Ra (b). The fitted curve is extracted using the experimental values. The numerical values are included for completeness.

FIG. 6.

Dependence of the non-dimensional critical time, Tc on Pe for various experiments with different Ra (a). Correlation of Tc with Pe and Ra (b). The fitted curve is extracted using the experimental values. The numerical values are included for completeness.

Close modal

In Fig. 6(b), we attempt to combine both the buoyancy contribution (Ra) and the advective-diffusive contribution (Pe) in an expression predicting the critical time, Tc, further extending the trends identified in the previous studies. The Tc values seem to collapse into the scaling: T c = P e 0.91 R a 0.55. The numerical data are not taken into account for the calculation of the relation. They appear to have slightly higher values than the experimental data although they qualitatively follow the same trend. This fact can be mainly attributed to the absence of any externally imposed disturbances in the numerical simulations, that indeed exist in experimental tries, such as vibrations or imperfect inlet conditions which may lead to the “premature” emergence of the instability in laboratory conditions. Nevertheless, certain perturbations in the simulations result from the numerical noise (e.g., due to the discrete grid). Depending on the physical parameters, the numerical noise may have a different influence at the delicate point of instability onset, causing the observed variation in the numerical data points.

Since a trend for the non-dimensional wavelength ( λ / h) of the instability could be identified from the shadowgraph images, another scaling is suggested in Fig. 7. The value of λ is defined as the average distance between two dark regions in the shadowgraph image across the cells' width (see embedded sketch in Fig. 7). In contrast with the previous studies that propose scalings using exclusively the advective flow contribution24 expressed by Pe or using exclusively the gravitational effect,36 here, it appears more appropriate to combine both contributing forces, resulting in the P e 0.23 R a 0.14 term. This reinforces the argument that two different mechanisms act to promote the onset of the instability. Ra is governed by Δ ρ, g, and h. The higher Ra is the more unstable the system gets and the wavelength decreases. On the other hand, a high Pe number means low diffusive mass transfer compared to the advective transfer by the injection flow. Thus, a higher driving density gradient is conserved, favoring the convective instability. In the opposite case of small Pe, the dynamics are strongly affected by diffusion, which assumes a stabilizing role (through erasing concentration gradients in small scales). Hence, the wavelength increases, or for very low Pe, the instability does not occur at all (cf. Fig. 5 for the critical Pe). The same tendency is observed for the numerical results ( in Fig. 7), with good qualitative and quantitative agreement. Slight deviations are due to the experimental conditions which can result in less pronounced density gradients, i.e., slight initial premixing at the inlet. This initial premixing further will vary to some extent with the experimental parameters, for example, the inlet velocity, contributing to the scatter of the experimental data points in Fig. 7.

FIG. 7.

Dependence of the non-dimensional wavelength, λ / h, on Pe and Ra. The fitted curve is extracted using the experimental values. The numerical values are included for completeness.

FIG. 7.

Dependence of the non-dimensional wavelength, λ / h, on Pe and Ra. The fitted curve is extracted using the experimental values. The numerical values are included for completeness.

Close modal

Compared to radial or laterally non-confined rectilinear assemblies, we also expect some discrepancies in the overall scalings of the occurring instability. The lateral walls induce a sheared (Poiseuille) flow and are also expected to confine the emerging flow structures. The manufacturing of the lateral walls, which consist of the PTFE spacers, might also slightly affect the flow and induce disturbances because of their nonzero roughness. In addition, the discrete h/w ratios used are expected to have a certain influence on the observed pattern wavelengths.

In some cases, a vortex exchange phenomenon is observed, termed as “tip splitting” in the previous literature16 for the radial case. This tip-splitting phenomenon was attributed to the radial divergence of the flow. Here, we observe a similar phenomenon of vortex generation and dissipation whose mechanism is not trivial to understand. As can be seen in Fig. 8, the counter-rotating vortices tend to diverge from their rectilinear direction and begin to position themselves toward the wall (Fig. 8, for t0). Figure 8(a) shows the experimental shadowgraph images for the case with Pe = 280 and Ra = 176 460. With regard to the “numerical shadowgraph” images, similar structures as in the experiments appear as shown in Fig. 8.

FIG. 8.

Progression of a displacement exhibiting the complex interaction between vortices for a case with Pe = 280 and Ra = 176 460 in experiment (a) and simulation (b). The red arrows signify the merge of two vortices.

FIG. 8.

Progression of a displacement exhibiting the complex interaction between vortices for a case with Pe = 280 and Ra = 176 460 in experiment (a) and simulation (b). The red arrows signify the merge of two vortices.

Close modal

The diverging rolls are probably due to the shear induced by the lateral walls and the non-uniform concentration gradients in the lateral direction, stemming from the Poiseuille profile itself. Later on, at t 0 + 2 s and t 0 + 4 s, the outermost vortices have already dissipated in Fig. 8(a) and the distance between the remaining vortices has increased. As a result of the diverging and dissipating rolls, there is space for the creation of a second generation of vortices in-between. The complex dynamics of the interaction between the first and second generation of vortices is difficult to decipher. In this case, there seems to be a fusion event toward the tail of the vortices, leaving only one vortex behind ( t 0 to t 0 + 2 s, noted in Fig. 8 with the red arrows). It is presumed that this vortex replacement will carry on as long as the density gradient across the displacement stratification is high enough to sustain the generation of new vortices replacing the dissipated ones. Diverging and newly forming rolls are likewise observed in the numerical results in Fig. 8(b). However, the vortex dissipation only sets in later ( t 0 + 14 s) due to the above-discussed delay in the pattern evolution. This kind of dynamic vortex interaction is not observed in previously studied rectilinear cases of the instability.24,32,36 We strongly suspect that the shear induced by the lateral walls ( d u x d y) might contribute to the phenomenology of the resulting vortex pattern. To further support this hypothesis, we simulated one case of the displacement using a free slip condition on the lateral walls. For this case, the vortex movement toward the sides disappears (for Pe = 280, Ra = 176 460, as described above), a fact that supports the hypothesis that the vortex–vortex–wall interactions appear because of the lateral shearing of the velocity profile. Note that a naturally diverging flow is present in the case of radial displacements studied experimentally16 causing vortex–vortex interactions of a different kind. In previous numerical studies in rectilinear geometries, where lateral walls were not present, no such vortex interactions were observed.40 

The above-described temporal pattern evolution is visually represented in Fig. 9 where the shadowgraph data from the forward part of the front are plotted against time (a moving sampling line is used, that progresses to the right with the average velocity of the displacement). The lateral movement of the vortices and the subsequent collision with the wall is visible, both for the experimental case [Fig. 9(a)] and the simulation with no-slip lateral walls [Fig. 9(b)]. The generation of new vortices can also be spotted for both cases in Figs. 9(a) and 9(b). For the simulation with free-slip lateral walls [Fig. 9(c)], the space–time plot shows that the divergent movement of the vortices is absent in this case, and no new vortices are generated. On the contrary, vortices disappear with time, leading to an increase in the pattern wavelength due to the gradually decreasing density gradient.

FIG. 9.

Space–time plot of the instability at the point behind the displacement's front tip for a case with Pe = 280 and Ra = 176 460 from shadowgraph experiment (a), numerical simulation (b), and numerical simulation with free-slip wall boundary conditions (c).

FIG. 9.

Space–time plot of the instability at the point behind the displacement's front tip for a case with Pe = 280 and Ra = 176 460 from shadowgraph experiment (a), numerical simulation (b), and numerical simulation with free-slip wall boundary conditions (c).

Close modal

To further quantify the instability, we investigated the flow field with μ-PIV. Only a small part of the experimental flow cell is captured in the μ-PIV recordings. In Fig. 10(b), the averaged flow field result is presented.

FIG. 10.

Schematic representation of the PIV measurements field of view marked by a white rectangle (a), and time-averaged velocity magnitude in the same area for the time the instability requires to pass through the field of view (approximately 25 s) (b).

FIG. 10.

Schematic representation of the PIV measurements field of view marked by a white rectangle (a), and time-averaged velocity magnitude in the same area for the time the instability requires to pass through the field of view (approximately 25 s) (b).

Close modal

The vector field is averaged throughout the time the instability is required to go through the entire optical field ( 25 s). In parallel, a shadowgraph image is shown, to better locate the flow field domain [Fig. 10(a)]. The flow field comprises the x and y velocity components. These values correspond to an averaged velocity magnitude over the z direction due to the depth of field resulting from the optical configuration; in this case, Δ z = 0.18 mm. This optical averaging only allows a rough representation of the velocity field. The two regions with high-velocity magnitude correspond to the two vortices that are close to the wall, in the region marked by a white rectangle in Fig. 10(a). This slight increase in the velocity magnitude is caused by an increase in the ux and uy component of the velocity vectors. The total movement of the eddies toward the wall might also contribute to the local increase of the uy component. On the contrary, the velocity magnitude decreases in the zones in-between the vortices because of the absence of velocity components in the y direction. The wall area is also visible in Fig. 10(b) for y < −0.5  mm where the velocity magnitude is close to zero.

To better visualize the progress of the velocity magnitude in time, Fig. 11(a) shows the temporal progression of the flow velocity magnitude at a stationary line for x 8 mm downstream of the flow cell entrance, in the sense of a space–time plot. In Fig. 11(b), a space–time plot for the numerical results representing the ux velocity magnitude at z = 0.2 mm, for x = 20 mm downstream of the inlet is shown.

FIG. 11.

Space–time plot for velocity magnitude for the case with Pe = 1850 and Ra = 123 930 (a): μ-PIV results representing velocity magnitude; (b) numerical results representing the ux, uy velocity magnitude for z = 0.2.

FIG. 11.

Space–time plot for velocity magnitude for the case with Pe = 1850 and Ra = 123 930 (a): μ-PIV results representing velocity magnitude; (b) numerical results representing the ux, uy velocity magnitude for z = 0.2.

Close modal

As is evident for the experimental case, the instability reaches the respective position in the flow cell at t  6  s. Before that point, there is no clear pattern visible in the flow velocity magnitude. The slight decrease in the velocity magnitude before the instability reaches x 8 mm is probably caused by a reduction in the seeding density of the PIV tracer particles due to some sedimentation before starting the injection of the second solution. After t  6 s, the existence of two instability vortices is clear. Toward the end, the existence of the vortices becomes less clear, as they start to dissipate, and as a result, the velocity magnitude becomes more uniform along the y–axis. In the current representation, we are observing the spatiotemporal progression of the instability. Specifically, we observe the progression of the counter-rotating vortices in time as we move toward the right-hand area in Fig. 11, but at the same time, the area of interest moves toward the upwind part of the vortices (tails) since we plot the data on a stationary line. In Fig. 11(b), the instability arises in the same way as in the experimental results. The magnitude of the velocities agrees well with the experimental results within the above-mentioned limitations of the PIV measurements. Two emerging vortices and their movement can be traced in the selected region.

In Fig. 12, the uz, uy velocity vector field and the magnitude for the uz (b) and uy (c) velocities are presented in a sectional view. A main couple of counterrotating vortices can be clearly seen in the center of the cross section of the flow geometry. Toward the two lateral walls, the vortices have a slightly different structure and the ascending components of the uz velocities appear more intense than in the center. This may compensate the descending gravity current in the advancing tip due to continuity. The outward movement of the outer vortices can also be observed as the vectors close to the bottom wall have a more pronounced vortex structure and higher velocity magnitudes pointing toward the lateral walls.

FIG. 12.

Numerical velocity vector field in a cross section of the flow domain 45 mm downstream the inlet for a case with Pe = 280 and Ra = 176 460 (a). The uz velocity magnitude is displayed by the background color in (b), while the uy velocity magnitude is displayed in (c).

FIG. 12.

Numerical velocity vector field in a cross section of the flow domain 45 mm downstream the inlet for a case with Pe = 280 and Ra = 176 460 (a). The uz velocity magnitude is displayed by the background color in (b), while the uy velocity magnitude is displayed in (c).

Close modal

In total, the cross section clearly illustrates the flow structure caused by the buoyancy instability of the lower side of the parabola-like profile (cf. Fig. 1), while the upper side is buoyantly stable and mainly contributes to the large-scale flow in the rectangular geometry.

In this study, we conducted for the first time a combined experimental and numerical study, providing qualitative and quantitative insights of the buoyancy-driven instability during the injection of a denser miscible fluid into a less dense one in a rectilinear laterally bounded geometry. We conducted experiments in a variety of gap widths, with different fluid density differences ( Δ ρ) and varied injection flow rates (q). We visualized the experiments using a shadowgraph technique and we conducted μ-PIV experiments for flow field calculations.

The dependence of the developing instability on the flow parameters, characterized by the dimensionless critical time, Tc, and the dimensionless wavelength, λ / h, was validated experimentally and numerically. Approximate scalings were suggested, combining the contribution of the advective flow via the Péclet number, and the buoyancy effect via the Rayleigh number. In more detail, it has been shown that the instability requires a minimum Pe value to occur, which decreases as the ScRa value increases. For small ScRa values, a stabilization of the system is observed. It has also been shown that as Pe and Ra increase the onset time of the instability, Tc decreases. The same effect is observed for the characteristic wavelength, λ / h. Complex vortex–vortex dynamics are presented for the rectilinear geometry, and a comparison with the free-slip case and radial cases from previous work is drawn. The effect of the lateral walls plays a significant role on the vortex dynamics. Initial PIV experiments provided the first quantitative information on the flow field. A comparison with the numerically obtained flow field supported the interpretation of the velocity measurements. The simulations provided further insights in the full three-dimensional (3D) velocity distribution which is not accessible in the experiments.

The extracted scalings of the studied instability can serve as a guide to predict or even suppress such instabilities in similar systems. This work forms a basis for future research on more general scaling laws, in the direction of a rigorous analysis, or extended governing parameters. The complex vortex–vortex interactions, e.g., the detailed dissipation mechanism, can further be elucidated as part of a future study.

See the supplementary material for more details on the flow cells geometry; the experimental parameters; the numerical simulation setup parameters (domain geometry, computational mesh); the methods to compute the critical time, Tc, and dimensionless wavelength λ / h; and the concentration distribution during the displacement resulting from the non-parabolic profile due to the gravity current.

Y. S. would like to thank Mert Uzun for his support with conducting the shadowgraph experiments. This work was supported by the German Aerospace Center (DLR) [Grant No. 50WM2061 (project ChemFront)]. A.T. and D.H. are grateful for the support of ESA (No. PRODEX 400138061) and National Research, Development, and Innovation Office (No. TKP2021-NVA-19). P.P. acknowledges the support by the New National Excellence Program of the Ministry for Culture and Innovation from the source of the National Research, Development, and Innovation Fund (No. ÚNKP-22-3-SZTE-441).

The authors have no conflicts to disclose.

Yorgos Stergiou: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Paszkal Papp: Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Dezso Horvath: Resources (equal); Supervision (equal); Writing – review & editing (equal). Agota Toth: Funding acquisition (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Kerstin Eckert: Resources (equal); Supervision (equal); Writing – review & editing (equal). Karin Schwarzenberger: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
H.
Erfani
,
N.
Karadimitriou
,
A.
Nissan
,
M. S.
Walczak
,
S.
An
,
B.
Berkowitz
, and
V.
Niasar
, “
Process-dependent solute transport in porous media
,”
Transp. Porous Media
140
,
421
435
(
2021
).
2.
C.
Cohen
,
M.
Berhanu
,
J.
Derr
, and
S.
Courrech du Pont
, “
Buoyancy-driven dissolution of inclined blocks: Erosion rate and pattern formation
,”
Phys. Rev. Fluids
5
,
053802
(
2020
).
3.
H. E.
Huppert
and
J. A.
Neufeld
, “
The fluid mechanics of carbon dioxide sequestration
,”
Annu. Rev. Fluid Mech.
46
,
255
272
(
2014
).
4.
N.
Heinemann
,
J.
Alcalde
,
J. M.
Miocic
,
S. J. T.
Hangx
,
J.
Kallmeyer
,
C.
Ostertag-Henning
,
A.
Hassanpouryouzband
,
E. M.
Thaysen
,
G. J.
Strobel
,
C.
Schmidt-Hattenberger
,
K.
Edlmann
,
M.
Wilkinson
,
M.
Bentham
,
R. S.
Haszeldine
,
R.
Carbonell
, and
A.
Rudloff
, “
Enabling large-scale hydrogen storage in porous media—The scientific challenges
,”
Energy Environ. Sci.
14
,
853
864
(
2021
).
5.
D.-W.
Chung
,
P. R.
Shearing
,
N. P.
Brandon
,
S. J.
Harris
, and
R. E.
García
, “
Particle size polydispersity in Li-ion batteries
,”
J. Electrochem. Soc.
161
,
A422
(
2014
).
6.
B.
Jha
,
L.
Cueto-Felgueroso
, and
R.
Juanes
, “
Fluid mixing from viscous fingering
,”
Phys. Rev. Lett.
106
,
194502
(
2011
).
7.
B.
Jha
,
L.
Cueto-Felgueroso
, and
R.
Juanes
, “
Synergetic fluid mixing from viscous fingering and alternating injection
,”
Phys. Rev. Lett.
111
,
144501
(
2013
).
8.
P.
Arosio
,
T.
Müller
,
L.
Mahadevan
, and
T. P. J.
Knowles
, “
Density-gradient-free microfluidic centrifugation for analytical and preparative separation of nanoparticles
,”
Nano Lett.
14
,
2365
2371
(
2014
).
9.
S.
Hassanpour Tamrin
,
A.
Sanati Nezhad
, and
A.
Sen
, “
Label-free isolation of exosomes using microfluidic technologies
,”
ACS Nano
15
,
17047
17079
(
2021
).
10.
J.
Nijhuis
,
S.
Schmidt
,
N. N.
Tran
, and
V.
Hessel
, “
Microfluidics and macrofluidics in space: ISS-proven fluidic transport and handling concepts
,”
Front. Space Technol.
2
, 779696 (
2022
).
11.
A.
Vailati
,
H.
Bataller
,
M. M.
Bou-Ali
,
M.
Carpineti
,
R.
Cerbino
,
F.
Croccolo
,
S. U.
Egelhaaf
,
F.
Giavazzi
,
C.
Giraudet
,
G.
Guevara-Carrion
,
D.
Horváth
,
W.
Köhler
,
A.
Mialdun
,
J.
Porter
,
K.
Schwarzenberger
,
V.
Shevtsova
, and
A.
De Wit
, “
Diffusion in liquid mixtures
,”
npj Microgravity
9
,
1
8
(
2023
).
12.
J. H. E.
Cartwright
,
B.
Escribano
,
C. I.
Sainz-Díaz
, and
L. S.
Stodieck
, “
Chemical-garden formation, morphology, and composition—II: Chemical gardens in microgravity
,”
Langmuir
27
,
3294
3300
(
2011
).
13.
S.
Chandrasekhar
,
Hydrodynamic and Hydromagnetic Stability
(
Courier Corporation
,
1981
).
14.
A.
De Wit
, “
Chemo-hydrodynamic patterns and instabilities
,”
Annu. Rev. Fluid Mech.
52
,
531
555
(
2020
).
15.
J.-R.
Authelin
,
F.
Brochard
, and
P.-G. D.
Gennes
, “ ‘
Instabilité gravitationnelle de deux fluides miscibles,‘ Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'univers
,”
Sci. de la Terre
317
,
1539
1541
(
1993
).
16.
F.
Haudin
,
L. A.
Riolfo
,
B.
Knaepen
,
G. M.
Homsy
, and
A.
De Wit
, “
Experimental study of a buoyancy-driven instability of a miscible horizontal displacement in a Hele-Shaw cell
,”
Phys. Fluids
26
,
044102
(
2014
).
17.
J.
Fernandez
,
P.
Kurowski
,
P.
Petitjeans
, and
E.
Meiburg
, “
Density-driven unstable flows of miscible fluids in a Hele-Shaw cell
,”
J. Fluid Mech.
451
,
239
260
(
2002
).
18.
F.
Graf
,
E.
Meiburg
, and
C.
Härtel
, “
Density-driven instabilities of miscible fluids in a Hele-Shaw cell: Linear stability analysis of the three-dimensional Stokes equations
,”
J. Fluid Mech.
451
,
261
282
(
2002
).
19.
P. G.
Saffman
and
G. I.
Taylor
, “
The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid
,”
Proc. R. Soc. London Ser. A
245
,
312
329
(
1958
).
20.
G. M.
Homsy
, “
Viscous fingering in porous media
,”
Annu. Rev. Fluid Mech.
19
,
271
311
(
1987
).
21.
R. J.
Glass
,
J.-Y.
Parlange
, and
T. S.
Steenhuis
, “
Wetting front instability—1: Theoretical discussion and dimensional analysis
,”
Water Resources Res.
25
,
1187
1194
(
1989
).
22.
R.
Krechetnikov
and
G. M.
Homsy
, “
On a new surfactant-driven fingering phenomenon in a Hele-Shaw cell
,”
J. Fluid Mech.
509
,
103
124
(
2004
).
23.
J.
Fernandez
,
R.
Krechetnikov
, and
G. M.
Homsy
, “
Experimental study of a surfactant-driven fingering phenomenon in a Hele-Shaw cell
,”
J. Fluid Mech.
527
,
197
216
(
2005
).
24.
A.
Mizev
,
E.
Mosheva
,
K.
Kostarev
,
V.
Demin
, and
E.
Popov
, “
Stability of solutal advective flow in a horizontal shallow layer
,”
Phys. Rev. Fluids
2
,
103903
(
2017
).
25.
L.
Rayleigh
, “
Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density
,”
Proc. London Math. Soc.
s1-14
,
170
177
(
1882
).
26.
G.
Taylor
, “
The instability of liquid surfaces when accelerated in a direction perpendicular to their planes—I
,”
Proc. R. Soc. London Ser. A
201
,
192
196
(
1950
).
27.
P.
Kurowski
,
C.
Misbah
, and
S.
Tchourkine
, “
Gravitational instability of a fictitious front during mixing of miscible fluids
,”
Eurphys. Lett.
29
,
309
(
1995
).
28.
P.
Garik
,
J.
Hetrick
,
B.
Orr
,
D.
Barkey
, and
E.
Ben-Jacob
, “
Interfacial cellular mixing and a conjecture on global deposit morphology
,”
Phys. Rev. Lett.
66
,
1606
1609
(
1991
).
29.
J. R.
Lister
and
R. C.
Kerr
, “
The propagation of two-dimensional and axisymmetric viscous gravity currents at a fluid interface
,”
J. Fluid Mech.
203
,
215
249
(
1989
).
30.
N.
Didden
and
T.
Maxworthy
, “
The viscous spreading of plane and axisymmetric gravity currents
,”
J. Fluid Mech.
121
,
27
42
(
1982
).
31.
P.
Petitjeans
and
T.
Maxworthy
, “
Miscible displacements in capillary tubes—Part 1: Experiments
,”
J. Fluid Mech.
326
,
37
56
(
1996
).
32.
S.
Obernauer
, “
Inestabilidades Entre Fluídos Miscibles en Medios Porosos
,” Ph.D. thesis (
Universidad de Buenos Aires. Facultad de Ciencias Exactas y Naturales
,
1999
).
33.
A.
Baker
,
Á.
Tóth
,
D.
Horváth
,
J.
Walkush
,
A. S.
Ali
,
W.
Morgan
,
Á.
Kukovecz
,
J. J.
Pantaleone
, and
J.
Maselko
, “
Precipitation pattern formation in the copper(II) oxalate system with gravity flow and axial symmetry
,”
J. Phys. Chem. A
113
,
8243
8248
(
2009
).
34.
Q.
Wang
,
K. S.
Hernesman
, and
O.
Steinbock
, “
Flow-driven precipitation patterns with microemulsions in a confined geometry
,”
ChemSystemsChem
2
,
e1900037
(
2020
).
35.
Y.
Stergiou
,
M. J. B.
Hauser
,
A.
De Wit
,
G.
Schuszter
,
D.
Horváth
,
K.
Eckert
, and
K.
Schwarzenberger
, “
Chemical flowers: Buoyancy-driven instabilities under modulated gravity during a parabolic flight
,”
Phys. Rev. Fluids
7
,
110503
(
2022
).
36.
G.
Pótári
,
Á.
Tóth
, and
D.
Horváth
, “
Precipitation patterns driven by gravity current
,”
Chaos
29
,
073117
(
2019
).
37.
E.
Balog
,
P.
Papp
,
Á.
Tóth
,
D.
Horváth
, and
G.
Schuszter
, “
The impact of reaction rate on the formation of flow-driven confined precipitate patterns
,”
Phys. Chem. Chem. Phys.
22
,
13390
13397
(
2020
).
38.
S.
Alqatari
,
T. E.
Videbaek
,
S. R.
Nagel
,
A. E.
Hosoi
, and
I.
Bischofberger
, “
Confinement-induced stabilization of the Rayleigh-Taylor instability and transition to the unconfined limit
,”
Sci. Adv.
6
,
eabd6605
(
2020
).
39.
L.
Talon
,
N.
Goyal
, and
E.
Meiburg
, “
Variable density and viscosity, miscible displacements in horizontal Hele-Shaw cells—Part 1: Linear stability analysis
,”
J. Fluid Mech.
721
,
268
294
(
2013
).
40.
M. O.
John
,
R. M.
Oliveira
,
F. H. C.
Heussler
, and
E.
Meiburg
, “
Variable density and viscosity, miscible displacements in horizontal Hele-Shaw cells—Part 2: Nonlinear simulations
,”
J. Fluid Mech.
721
,
295
323
(
2013
).
41.
K. S.
Gage
and
W. H.
Reid
, “
The stability of thermally stratified plane Poiseuille flow
,”
J. Fluid Mech.
33
,
21
32
(
1968
).
42.
A. A.
Mohamad
and
R.
Viskanta
, “
Laminar flow and heat transfer in Rayleigh–Bénard convection with shear
,”
Phys. Fluids A
4
,
2131
2140
(
1992
).
43.
N.
Goyal
and
E.
Meiburg
, “
Miscible displacements in Hele-Shaw cells: Two-dimensional base states and their linear stability
,”
J. Fluid Mech.
558
,
329
355
(
2006
).
44.
R. M.
Oliveira
and
E.
Meiburg
, “
Saffman-Taylor instability and the inner splitting mechanism
,”
Phys. Rev. Lett.
118
,
124502
(
2017
).
45.
W.
Merzkirch
,
Flow Visualization
(
Academic Press
,
1987
).
46.
M.
Berhanu
,
J.
Philippi
,
S.
Courrech du Pont
, and
J.
Derr
, “
Solutal convection instability caused by dissolution
,”
Phys. Fluids
33
,
076604
(
2021
).
47.
J.
Philippi
,
M.
Berhanu
,
J.
Derr
, and
S.
Courrech du Pont
, “
Solutal convection induced by dissolution
,”
Phys. Rev. Fluids
4
,
103801
(
2019
).
Published open access through an agreement with Technische Universitat Dresden Institut fur Verfahrenstechnik und Umwelttechnik

Supplementary Material