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, *T _{C}*, and dimensionless wavelength, $ \lambda / 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.

## I. INTRODUCTION

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} CO_{2} sequestration,^{3} H_{2} 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 occur^{13–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 $ \mu 1 \u2265 \mu 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 instability^{27} as depicted in Fig. 1.

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 Kerr^{29} and Didden and Maxworthy^{30} 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. Obernauer^{32} 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 particles^{33} in the contact zone of two reactant solutions during radial displacement. Patterns in confined flow-driven microemulsion generation^{34} 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 studies^{16,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).

## II. PHYSICAL FORMULATION AND GOVERNING EQUATIONS

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 $ \rho 2 > \rho 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}

^{43,44}the phenomena are governed by the following three-dimensional, incompressible Navier–Stokes equations:

*ρ*is the density,

*μ*is the viscosity, and $ g \u2192$ is the gravitational acceleration.

*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.

Let us then define the characteristic scales as $ L * = h , \u2009 T * = h 2 / D , \u2009 U * = L * / T *$, and $ P * = D \mu / h 2$, and the dimensionless quantities *ζ*, $ t \u0303 , \u2009 u \u0303 , \u2009 p \u0303$.

*Sc*and

*Ra*.

## III. METHODS

### A. Experimental setup

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 experiments^{26,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.

Fluid . | Glycerol (% v/v) . | Density ( $ g / cm 3$) . | Viscosity ( $ cP$) . |
---|---|---|---|

1 | 0 | 0.9972 | 0.92 |

2 | 6.8 | 1.0182 | 1.2209 |

3 | 13.6 | 1.0363 | 1.5445 |

Fluid . | Glycerol (% v/v) . | Density ( $ g / cm 3$) . | Viscosity ( $ cP$) . |
---|---|---|---|

1 | 0 | 0.9972 | 0.92 |

2 | 6.8 | 1.0182 | 1.2209 |

3 | 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 \u2009 \mu 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 $ \u2248 0.18 \u2009 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 \u2009 kHz$.

### B. Numerical simulations

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 \u2212 9 \u2009 m 2 \u2009 s \u2212 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 $ \Delta t = 0.5 \u2026 2 \u2009 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 $ \Delta x = 0.2 \u2009 mm , \u2009 \Delta y = 0.05 \u2009 mm , and \u2009 \Delta z = 0.04 \u2009 mm$. A grid convergence study was performed with a refined mesh (50% higher resolution for $ \Delta x , \u2009 \Delta y$, and $ \Delta z$, and a time step, $ \Delta 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 $ \u2207 2 c z$ field is shown, with *c _{z}* being the

*z*-averaged concentration values.

## IV. RESULTS AND DISCUSSION

### A. Shadowgraph observations: Characteristics of the instability

In Fig. 3, shadowgraph experiments of a $ 6.8 \u2009 %$ glycerol solution injected into distilled water in a flow cell ( $ h = 0.8 \u2009 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.

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 $ \Delta \rho $ between the two fluids, the influence of density gradient is apparent (Fig. 4). Larger $ \Delta \rho $ 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.

A similar trend for the instability wavelength depending on *h* has been previously observed for experiments with radial injection^{16} 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 $ \Delta \rho $ 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}

### B. Scaling of the instability onset

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 $ \Delta \rho $, *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 = \Delta \rho g h 3 / \rho D 2$. In the chart, $\xb0$ 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 $ \Delta \rho $) 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 \u221d ( ScRa ) \u2212 0.39$. We can, thus, assume that the critical condition in our experiments is mainly controlled by the characteristic length (*h*), the density difference ( $ \Delta \rho $), and the magnitude of the advective flow (*U*). All of the above can be combined using the quantities expressed by *ScRa* and *Pe*.

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, *t _{c}*, 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,

*t*, the standard deviation of the grayscale value of the shadowgraph images (or the $ \u2207 2 c z$, for the CFD results) is used. The

_{c}*t*is approximately set to the time when the standard deviation starts increasing rapidly (cf. see supplementary material).

_{c}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 $ \Delta \rho $, while the only varying parameter affecting *Pe* was the flow rate, *q. T _{c}* shows a weak proportional decrease according to $ \u223c 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 instability

^{38}or for a solutal system, based on the Rayleigh–Bénard scaling laws.

^{46,47}Although other studies suggest an independency of

*t*by

_{c}*Pe*,

^{16,36}we suppose this is caused by the relatively high flow rates used in those cases, where the differences between

*t*might be negligible.

_{c}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, *T _{c}*, further extending the trends identified in the previous studies. The

*T*values seem to collapse into the scaling: $ T c = P e \u2212 0.91 R a \u2212 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.

_{c}Since a trend for the non-dimensional wavelength ( $ \lambda / 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 contribution^{24} 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 \u2212 0.23 R a \u2212 0.14$ term. This reinforces the argument that two different mechanisms act to promote the onset of the instability. *Ra* is governed by $ \Delta \rho $, *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 ( $\u2666$ 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.

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.

### C. Pattern evolution

In some cases, a vortex exchange phenomenon is observed, termed as “tip splitting” in the previous literature^{16} 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 *t*_{0}). 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.

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 \u2009 s$ and $ t 0 + 4 \u2009 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 \u2009 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 \u2009 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 experimentally^{16} 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.

### D. Velocity field

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.

The vector field is averaged throughout the time the instability is required to go through the entire optical field ( $ \u2248 25 \u2009 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, $ \Delta z = 0.18 \u2009 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 *u _{x}* and

*u*component of the velocity vectors. The total movement of the eddies toward the wall might also contribute to the local increase of the

_{y}*u*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}*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 \u2248 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 *u _{x}* velocity magnitude at

*z*= 0.2 mm, for

*x*= 20 mm downstream of the inlet is shown.

As is evident for the experimental case, the instability reaches the respective position in the flow cell at $ t \u2009 \u2248$ 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 \u2248 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 \u2248$ 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 *u _{z}*,

*u*velocity vector field and the magnitude for the

_{y}*u*(b) and

_{z}*u*(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

_{y}*u*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.

_{z}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.

## V. CONCLUDING REMARKS

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 ( $ \Delta \rho $) 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, *T _{c}*, and the dimensionless wavelength, $ \lambda / 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,

*T*decreases. The same effect is observed for the characteristic wavelength, $ \lambda / 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.

_{c}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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Hydrodynamic and Hydromagnetic Stability*