Experimental and numerical investigation of a density-driven instability during a horizontal miscible displacement

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,¹ geochemical flows,² CO₂ sequestration,³ H₂ storage,⁴ or Li-ion battery design.⁵ 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.¹² 

It has been observed that when a fluid 1 of density ρ₁ and viscosity μ₁ is injected into a miscible fluid 2 of density ρ₂ and viscosity μ₂, 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,²¹ 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 instability²⁷ as depicted in Fig. 1.

Garik et al.²⁸ 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²⁹ and Didden and Maxworthy³⁰ 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,³¹ but due to spatial confinement, the phenomenon did not produce more than one set of instability rolls in the flow cell. Obernauer³² 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³³ in the contact zone of two reactant solutions during radial displacement. Patterns in confined flow-driven microemulsion generation³⁴ have been likewise attributed to such convective instabilities. Later on, Haudin et al.¹⁶ 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.³⁵ In continuation of this work,¹⁶ Pótári et al.³⁶ described the phenomenon analytically and coupled it to experimental observations present in precipitation patterns. The same system has been revisited by Balog et al.³⁷ Subsequently, scaling laws for high Pe number cases have been proposed.³⁶ 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.³⁸ Numerical simulations by Talon et al.³⁹ and John et al.⁴⁰ 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).

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 ρ₁. 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,¹⁶ 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.⁴² 

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

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).⁴⁵ 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 c_z 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.

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.

A similar trend for the instability wavelength depending on h has been previously observed for experiments with radial injection¹⁶ and numerical simulations in a rectilinear geometry with free slip walls.³⁶ 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.,⁴⁰ 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.³⁶ 

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.

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.³⁸ 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_c, the standard deviation of the grayscale value of the shadowgraph images (or the $∇ 2 c z$⁠, for the CFD results) is used. The t_c 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. T_c shows a weak proportional decrease according to $∼ 1 / P e$⁠, a finding in accordance with previous studies in similar systems.²⁴ For systems without advection, scaling laws of similar form have been proposed: either for the stationary Rayleigh–Taylor instability³⁸ or for a solutal system, based on the Rayleigh–Bénard scaling laws.^46,47 Although other studies suggest an independency of t_c by Pe,^16,36 we suppose this is caused by the relatively high flow rates used in those cases, where the differences between t_c might be negligible.

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_c 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 contribution²⁴ expressed by Pe or using exclusively the gravitational effect,³⁶ 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.

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 literature¹⁶ 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₀). 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 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 experimentally¹⁶ 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.⁴⁰ 

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.

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 (⁠ $≈ 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 u_x and u_y component of the velocity vectors. The total movement of the eddies toward the wall might also contribute to the local increase of the u_y 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 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 ≈$ 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 u_z, u_y velocity vector field and the magnitude for the u_z (b) and u_y (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 u_z 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.

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, T_c, 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, T_c 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.