Concentration-based velocity reconstruction in convective Hele–Shaw flows

Abstract

We examine the process of convective dissolution in a Hele–Shaw cell. We consider a one-sided configuration and we propose a newly developed method to reconstruct the velocity field from concentration measurements. The great advantage of this Concentration-based Velocity Reconstruction (CVR) method consists of providing both concentration and velocity fields with a single snapshot of the experiment recorded in high resolution. We benchmark our method vis–à–vis against numerical simulations in the instance of Darcy flows, and we also include dispersive effects to the reconstruction process of non-Darcy flows. The absence of laser sources and the presence of one low-speed camera make this method a safe, accurate, and cost-effective alternative to classical PIV/PTV velocimetry processes. Finally, as an example of possible application, we employ the CVR method to analyse the tip splitting phenomena.

Graphic abstract

1 Introduction

The Hele–Shaw cell consists of two paralleled and transparent plates located within a uniform gap. When the gap is sufficiently thin, the flow inside the cell can reproduce a Stokes flow. This feature together with optical accessibility represents great advantages of the Hele–Shaw cell, making it a suitable tool to investigate in detail different flow configurations (De Wit 2016). In this work, we will focus on buoyancy-driven Hele–Shaw flows. The driving force of the flow consists of density differences induced by the presence of a solute and, more precisely, by the solute concentration gradients existing in the domain. In these conditions, under certain assumptions, the Hele–Shaw apparatus may accurately mimic a Darcy flow in a porous medium. This problem is of practical interest for a number of geophysical and environmental applications such as water contamination (Molen and Ommen 1988; LeBlanc 1984), carbon sequestration (Huppert and Neufeld 2014; Emami-Meybodi et al. 2015), sea ice formation (Feltham et al. 2006; Wettlaufer et al. 1997; Middleton et al. 2016), and granular flows (Lange et al. 1998), to name a few.

Velocity measurements in these flows are of crucial importance for the characterisation of the flow patterns, and it is a very active topic of research (Ehyaei and Kiger 2014; Kreyenberg et al. 2019; Kislaya et al. 2020; Anders et al. 2020). While numerical simulations can give detailed information in terms of both velocity and concentration fields, Hele–Shaw experiments provide qualitative (Fernandez et al. 2002; Salibindla et al. 2018; Thomas et al. 2018; Lemaigre et al. 2013; Pringle et al. 2002) and quantitative observations of the concentration field (Slim et al. 2013; De Paoli et al. 2020). However, accurate experimental measurements of the Eulerian velocity field are hard to obtain. In this work, we aim precisely at bridging this gap. We propose a new method to simultaneously reconstruct the concentration and velocity distributions in buoyancy-driven Hele–Shaw flows. We investigate experimentally the evolution of a buoyancy-driven flow in a one-sided configuration (Hewitt et al. 2013; De Paoli et al. 2017). We consider a rectangular domain in which the fluid density is maximum at the top wall, so that the configuration is unstable. Density differences existing across the fluid layer are induced by the presence of a solute that, dissolving and mixing in time, controls the evolution of the flow.

We developed a Concentration-based Velocity Reconstruction (CVR) method to accurately estimate the Eulerian velocity field in Hele-Shaw flows. The method is solely based on the fields obtained from the concentration measurements and relies on two main steps. First, we record the experiment by taking images of the solute distribution and we infer the solute concentration field over the cell. Then, we use the concentration distribution to solve the momentum equation of the flow and we reconstruct the Eulerian velocity field. The flow inside the Hele–Shaw cell is usually considered as an ideal Darcy flow. However, when the strength of viscous and diffusive contributions is comparable to that of convection, additional non-Darcy effects come into play and may considerably change the evolution of the flow. As a result, corrections to the ideal Darcy flow should be taken into account (Letelier et al. 2019; De Paoli et al. 2020). We show that the CVR algorithm can be easily adapted to non-Darcy flows and we employ this method to analyse the tip splitting phenomenon.

The paper is organised as follows. In Sect. 2, we describe the experimental setup. In Sect. 3, we introduce the CVR algorithm and we benchmark our method vis-à-vis against numerical Darcy simulations. Finally, in Sect. 4, we employ the CVR to analyse the problem of tip splitting.

2 Experimental setup

We consider the process of convective dissolution in a Hele–Shaw cell. In particular, we chose a one-sided configuration, consisting of a rectangular domain initially filled with pure fluid. All boundaries are impenetrable to fluid and solute, but the top boundary, where the solute concentration is kept constant during the experiments. As a result, the solute dissolves and, in correspondence of the top boundary, convective structures form and control the evolution of the flow. The system is sketched in Fig. 1. We record the evolution of the flow and we infer density fields at high resolution from transmitted light intensity (Slim et al. 2013; Ching et al. 2017; De Paoli et al. 2020).

The cell consists of two glass plates $(10\times 140\times 250 \,\hbox {mm}^{3})$ separated by rubber seals of different thickness, b, which define the dimensions of the fluid layer (width L and height H). The material used for the sealings is a high-quality impermeable rubber, based on aramid fibre with nitrile binder (Klinger–sil C–4400). Metal shims are located between the plates, which are held in place by an external metal framework and a set of 16 screws. The torque applied to all screws is constant to ensure the uniformity of the gap. Along the top wall, lies a steel mesh ($100\,\mu \hbox {m}$ grid size) that contains the solute powder. The cell is illuminated from the backside by a tunable system, consisting of an array of 150 LED lamps covered by an opaque glass, which makes the light distribution uniform over the cell. In the following, we will describe the fluids adopted and the processes of concentration reconstruction.

Fig. 1

Experimental setup. The Hele–Shaw cell is reported, with explicit indication of domain dimensions (L, H), reference frame (x, z), and main components of the apparatus, consisting of camera, pump, illumination system. Along the steel mesh located at the top wall, the dye is poured, and therefore, the mixture is considered as saturated (i.e., solute concentration is constant, $C=C_{s}$)

2.1 Working fluids

An aqueous solution of potassium permanganate (KMnO4) and water are used to mimic, respectively, the denser and lighter fluids. In the present configuration, the cell is initially filled with water (pure fluid), i.e.,, the solute concentration C at time $t=0$ is $C(x,z,t=0)=0$, where x, z are the spatial coordinates in horizontal and vertical directions, respectively. The dissolution of KMnO4 takes place just below the grid on top of the cell (i.e., at $z=H$, with H fluid layer height). Therefore, we assume that at the top boundary, the fluid is saturated of solute and the concentration, C, is constant during the experiment [$C(x,z=H,t)=C_{s}$, being $C_{s}$ the saturation value]. Viscosity and diffusion coefficient are assumed constant and independent of the solute concentration (Slim et al. 2013), and are, respectively, $\mu =9.2\times 10^{-4}\,\hbox {Pa}\cdot \hbox {s}$ and $D=1.65\times 10^{-9}\,\hbox {m}^{2}/\hbox {s}$. Based on the correlations proposed by Novotný and Söhnel (1988), we assume that the density of an aqueous solution of KMnO4 depends on fluid temperature, $\vartheta$, water density, $\rho (0)$, and KMnO4 concentration. In particular, we measure the fluid temperature and we use it to estimate the water density. The density difference between the saturated solution and pure water is, $\rho _{s}-\rho (0)=38\,\hbox {kg}/\hbox {m}^{3}$, being $\rho _{s}=\rho (C=C_{s})$ and $C_{s}=46\,\hbox {kg}/\hbox {m}^{3}$ the effective saturation value of concentration. The density of the mixture, $\rho (C)$, can be well approximated as a linear function of the solute concentration, as shown in Fig. 2a, by the equation:

$$\begin{aligned} \rho =\rho _s\biggl [1+\frac{\Delta \rho }{\rho _sC_s}\bigl (C-C_s\bigr )\biggr ]\textit{ ,} \end{aligned}$$

(1)

with $\Delta \rho =\rho _{s}-\rho (0)$.

Grains of potassium permanganate (grain size greater than $200\, \mu \hbox {m}$) are poured on the grid sitting on top of the cell and are pressed to form a compact layer. In this way, the porosity of the layer is reduced and, as a result, it will be difficult for the liquid subsequently injected to form a liquid layer on top of the grid. The light fluid (water) is injected with the aid of a syringe pump from two channels located at the bottom of the cell and the fluid level is increased up to the upper boundary. We remark here that it is crucial that the level of water has to reach the height at which the grid is located, and the fluid has to get in contact with it along the whole cell width. During the first contact, water is sucked by the solute grains, which are initially dry. Afterwards, the water content of the solid powder increases, keeping the KMnO4 layer solid but humid. Finally, when the powder layer is saturated with water, the suction process stops and the dissolution of KMnO4 in water takes place. After the conclusion of this phase, which takes about 2 s from first contact to powder saturation, the pump is shut down. Solute dissolves in water and a fluid layer denser than water forms below the grid: After an initial diffusive phase, the layer of heavy mixture thickens, becomes unstable, and finger like structures form (Slim 2014). The onset of convection may start earlier than expected (i.e., fingers form earlier) in the presence of local perturbations of solute concentration (Riaz et al. 2006; Riaz and Cinar 2014, and references therein) and interface movement (Myint and Firoozabadi 2013). We preformed experiments in order to limit the solute and interface perturbations induced by the shutdown of the pump, which have an effect only on the initial flow evolution (Slim 2014).

Fig. 2

a Variation of the mixture density, $\rho (C)$, with respect to pure water, $\rho (0)$, as a function of the solute concentration, C. The density value obtained with the polynomial correlation by Novotný and Söhnel (1988) (symbols) is very well fitted by a linear function [Eq. (1), solid line]. The density profile is shown here in correspondence of a temperature $\vartheta =25^{\circ } C$. b Mass fraction—light intensity calibration curve for a cell gap $b=0.30~\hbox {mm}$ and for a constant illumination condition (i.e., constant tension applied to the LED lamps). The solute mass fraction $\omega$ is shown as a function of the normalised light intensity $I(\omega )/I({0})$, i.e., the light intensity measured for a mixture with mass fraction $\omega$ divided by the light intensity measured for pure water. The maximum normalised light intensity value, $I(\omega )/I({0})=1$, corresponds to the minimum solute mass fraction ($\omega =0$, pure water). Experimental measurements (symbols) are well fitted by a third-order polynomial (solid line). Uncertainties on the light intensity and mass fraction are also shown (details are provided in Sect. 3.3)

2.2 Reconstruction of the concentration field

The solute concentration field is inferred from light intensity distribution. A flowchart of the reconstruction process is shown in Fig. 3. First, the fluid temperature is measured and it is used to compute water density, $\rho (0)$. The evolution of the flow is recorded by the camera (Canon EOS 1300D, $3456 \times 5184\,\hbox {px}$). The collected images, which contain the value of light intensity in each pixel, $I_{0}(px,px)$, are then corrected with basic image processing operations (conversion to grey scale, correction of trapezoidal distortion) to obtain the corrected light distribution, I(px, px). Finally, the spatial calibration is applied and the corrected light intensity distribution over the cell, I(x, z), is found.

At this stage, the solute mass fraction field is determined using the intensity-mass fraction calibration data. The calibration process is summarised as follows. We prepared a number of samples of KMnO4 aqueous solution, corresponding to different values of the solute mass fraction $\omega$, from pure water to saturation. To this aim, we used a high-precision scale (Sartorius Acculab Atilon model ATL-423-I, $\pm 10^{-4}\,\hbox {g}$). The cell is then filled with each fluid sample and the light intensity distribution is measured by the camera. In correspondence of each sample, the mean light intensity is computed. The calibration curve obtained in correspondence of fixed tension applied to the LEDs and gap width is shown in Fig. 2b, where the solute mass fraction $\omega$ is reported as a function of the normalised light intensity, i.e., the space-average light intensity measured for the sample, $I(\omega )$, divided by the space-averaged light intensity measured for pure water, I(0) (we defined as space average as the average taken over the fluid layer). We remark here that a proper calibration is required for each value of the cell gap and illumination condition (Alipour and De Paoli 2019). With the aid of the calibration curves $\omega (I)$ (fitted here via third-order polynomials), the solute mass fraction distribution, $\omega (x,z)$, is reconstructed over the entire image. Since the solute concentration is defined as $C(x,z)=\rho (C)\omega (x,z)$ and given Eq. (1), the concentration distribution over the cell is finally inferred.

The concentration field C(x, z) may be subsequently used as input for the Concentration-based Velocity Reconstruction algorithm, as well as to estimate other quantities of the flow, such as the solute dissolution rate (De Paoli et al. 2020), the mean scalar dissipation rate (De Paoli et al. 2019; Hidalgo et al. 2012), or the mixing length (De Wit 2004; Gopalakrishnan et al. 2017).

Fig. 3

Process of Concentration-based Velocity Reconstruction. The algorithm consists of two parts: the concentration reconstruction (yellow box) and the velocity reconstruction (grey box). First, the light intensity distribution, I(px, px), is taken from the camera at a frame rate of 1fps. Using spatial and intensity calibration data, the spatial distribution of the solute mass fraction, $\omega (x,z)$, is obtained. The fluid temperature $(\vartheta )$ is used to find the correct pure water density, $\rho (0)$, and the concentration distribution C(x, z) is obtained with the correlations proposed by Novotný and Söhnel (1988). The solute concentration distribution is then used as input for the velocity reconstruction: after finding the stream function, $\psi (x,z)$, from the momentum equation, the velocity components, u and w, are determined. The procedure is iterated for all snapshots considered (time loop) to find the time-dependent evolution of the flow. However, one frame is sufficient to reconstruct the instantaneous velocity field

3 Concentration-based velocity reconstruction (CVR)

Consider the system sketched in Fig. 4, where a side view of the Hele–Shaw cell is proposed. We refer to our experimental setup, in which the fluid properties ($\Delta \rho ,\mu ,D$) are fixed. The only property that can be varied is the gap thickness, b. When the gap thickness is small, ideally infinitesimal, the flow in the Hele–Shaw cell is a Stokes flow and reproduces a two-dimensional Darcy flow (Fernandez et al. 2002; Slim et al. 2013; Oltean et al. 2008). More precisely, this condition occurs when the Reynolds number of the flow (${\textit{Re}}$, based on the gap-averaged velocity and cell thickness) is ${\textit{Re}}\rightarrow 0$. However, when a solute is present, density-driven convection represents the driving force of the flow and the system is controlled by the complex interplay of vertical advection and transverse (i.e., wall-normal) diffusion. In particular, if transverse diffusion dominates over vertical advection, the velocity at which the solute diffuses in y direction is much higher than the vertical advective velocity. The consequent wall-normal concentration profile is nearly flat (Fig. 4a), and the system can still be considered two-dimensional and controlled by the Darcy law. We define this regime as Darcy flow or Darcy regime. If the cell gap is increased, the hydrodynamic resistance experienced by the flow decreases, with a consequent increase of the velocity. As a result, the wall-normal velocity gradient may induce an additional solute dispersion, also defined as Taylor or hydrodynamic dispersion (Taylor 1953). In this situation, a concentration gradient exists across the cell and the flow field is not anymore well described by a Darcy model. However, the flow can still be considered two-dimensional provided that the effect of solute dispersion induced by the presence of the walls is taken into account (Fig. 4b). In this case, the flow is defined as a Hele–Shaw flow. Letelier et al. (2019) have shown theoretically and numerically that corrections can be applied to the Darcy equation to recover the additional solute spreading induced by the hydrodynamic dispersion. This has been also confirmed experimentally in a recent study by De Paoli et al. (2020). Finally, if the cell gap is further increased, the formation of a second finger takes place (Fig. 4c), and the system has a full three-dimensional character (three-dimensional regime).

To identify the flow regime, a quantitative estimation of the relative importance of vertical advection, transverse diffusion, and hydrodynamic dispersion is required, and two-dimensionless parameters are considered: the Rayleigh number and the anisotropy ratio, respectively:

$$\begin{aligned} {\textit{Ra}}=\frac{g\Delta \rho b^{2}H}{12\mu D}\qquad \text { and }\qquad \epsilon =\frac{b}{\sqrt{12}H} . \end{aligned}$$

(2)

The first, ${\textit{Ra}}$, measures the relative strength of vertical advection and diffusion along the cell height. The second governing parameter, $\epsilon$, scales as the characteristic time of transverse diffusion $(b^{2}/D)$ and longitudinal advection $(H/\overline{w})$, where $\overline{w}$ is the gap-averaged vertical velocity. The combination of these two quantities gives a further dimensionless parameter, $\epsilon ^{2}{\textit{Ra}}$, which determines the flow behaviour (Letelier et al. 2019). When the gap thickness is infinitesimal ($b\rightarrow 0$), we have that $\epsilon ^{2}{\textit{Ra}}\rightarrow 0$: diffusion in wall-normal direction is faster then longitudinal convection, and the flow is purely two-dimensional and well described by the Darcy model (Fig. 4a). When b in increased, the flow is influenced by the presence of the walls, which introduce an additional solute dispersion. However, the behaviour of the flow can still be considered two-dimensional, and only one finger is present in y direction (Fig. 4b). In this case, we have that $\epsilon ^{2}{\textit{Ra}}\ll$. Finally, if b is further increased and $\epsilon ^{2}{\textit{Ra}}\ge 1$, the flow assumes a three-dimensional character (Fig. 4c): more than one finger forms in wall-normal direction and the system cannot be considered uniform across the gap anymore. In Sect. 3.1, we describe the method of velocity reconstruction in a Hele–Shaw cell in the presence of a Darcy flow, i.e., when $\epsilon ^{2}{\textit{Ra}}\rightarrow 0$. When $\epsilon ^{2}{\textit{Ra}}\ll 1$ (Hele–Shaw flow), the effect of dispersion has to be taken into account via corrective terms to the Darcy equation. However, the Concentration-based Velocity Reconstruction can still be applied, and it is presented in detail in Appendix A.

Fig. 4

Sketch of solute distribution across the cell gap in the three regimes: (a) Darcy regime, $\epsilon ^{2}{\textit{Ra}}\rightarrow 0$, (b) Hele–Shaw regime, $\epsilon ^{2}{\textit{Ra}}\ll 1$, and (c) three-dimensional regime, $\epsilon ^{2}{\textit{Ra}}\ge 1$

3.1 Velocity reconstruction of Darcy flows $\left( \epsilon ^{2}{\textit{Ra}}\rightarrow 0\right)$

We are in the presence of a Darcy flow if transverse diffusion dominates over convection (i.e., $\epsilon ^{2}{\textit{Ra}}\rightarrow 0$), and the momentum equation is described by the Darcy law. We consider a Newtonian, incompressible fluid and we assume that the Boussinesq approximation can be applied.¹ Therefore, the two-dimensional velocity field is controlled by the equations:

$$\begin{aligned} \nabla \cdot \mathbf {u}= \,& {} 0 \end{aligned}$$

(3)

$$\begin{aligned} \mathbf {u}= & {} -\frac{b^{2}}{12\mu }(\nabla p+\rho \mathbf {g})\text { ,} \end{aligned}$$

(4)

where $\mathbf {u}$ is the Darcy velocity (i.e., the gap-average velocity), p is the pressure, and $\mathbf {g}$ is the acceleration due to gravity (see electronic supplemental material for the derivation of the Darcy equation from the Navier–Stokes solution of a Poiseuille flow). We consider the fluid density as the only physical property that depends on concentration through the Eq. (1), and from Eq. (4), we end up with:

$$\begin{aligned} \mathbf {u}=-\frac{b^{2}}{12\mu }\left[ \nabla P+\frac{\Delta \rho }{C_{s}}C\mathbf {g}\right] \text { ,} \end{aligned}$$

(5)

where $P=p-\rho (0)gz$ is the reduced pressure. The procedure of velocity reconstruction is based on the information available from the experimental measurements of concentration, assuming that the system is controlled by Eqs. (3) and (5).

We introduce the stream function $\psi$ (Batchelor 1968, p. 78) defined as $(u,w)=(\partial \psi /\partial z,-\partial \psi /\partial x)$. Taking the curl of Eq. (5) and assuming two-dimensional flow, we correlate the concentration gradient to the stream function via the Poisson equation:

$$\begin{aligned} \nabla ^{2} \psi =\alpha \frac{\partial C}{\partial x}, \end{aligned}$$

(6)

being $\alpha =b^{2}g\Delta \rho /(12\mu C_{s})$. The right-hand side of Eq. (6) is obtained from the experimental measurements and the problem is now formulated in terms of stream function, $\psi (x,z)$.

Fig. 5

Boundary conditions for concentration (a) and velocity and stream function (b). In particular, $\psi =0$ on the boundaries is used to solve Eq. (6), based on the experimental configuration adopted. The reference frame (x, z) is also indicated, while $\mathbf {n}$ stands for the unit vector normal to the boundary

The system is sketched in Fig. 5, with indication of the boundary conditions used for concentration (5a), velocity, and stream function (5b). At the top boundary, the concentration is constant $[C(x,z=H)=Cs]$ and no-penetration condition is applied in z-direction [$w(x,z=H)=0$]. Equation (6) and definition of $\psi$ give:

$$\begin{aligned} \nabla ^{2} \psi (x,z=H)=0\qquad \text { and} \qquad \frac{\partial \psi }{\partial x}\bigg |_{(x,z=H)}=0\text { .} \end{aligned}$$

(7)

We assume no flux of solute ($\partial C/\partial x=0$) and no penetration ($u=0$) across the side walls, which give:

$$\begin{aligned} \nabla ^{2} \psi (x=0,z)= & {} 0\qquad \text { and} \qquad \frac{\partial \psi }{\partial z}\bigg |_{(x=0,z)}=0\text { ,} \end{aligned}$$

(8)

$$\begin{aligned} \nabla ^{2} \psi (x=L,z)= & {} 0\qquad \text { and} \qquad \frac{\partial \psi }{\partial z}\bigg |_{(x=L,z)}=0\text { } \end{aligned}$$

(9)

for the left and right boundaries, respectively. Finally, at the bottom boundary, no-penetration condition ($w=0$) gives:

$$\begin{aligned} \frac{\partial \psi }{\partial x}\bigg |_{(x,z=0)}=0. \end{aligned}$$

(10)

From Eqs. (7)–(10), we conclude that $\psi$ is constant along the boundaries, and since the streamlines are impenetrable to fluid, all boundaries belong to the same streamline. We arbitrarily set the value of the stream function on the boundaries to zero, so that

$$\begin{aligned} \psi =0 \text { on }\Omega \text { ,} \end{aligned}$$

(11)

being $\Omega$ the boundary of the system. The Concentration-based Velocity Reconstruction (CVR) proposed here consists of solving numerically Eq. (6), where the r.h.s. is determined experimentally and the boundary conditions applied are obtained from (7)–(11). To this aim, we employ a sixth-order nine-point compact finite-difference scheme (Nabavi et al. 2007), where velocity components and all the other spatial derivatives are estimated using a sixth-order finite-difference scheme.

3.2 Algorithm validation with numerical results (Darcy flow)

The CVR algorithm is initially validated with data obtained from numerical, two-dimensional Darcy simulations. To this aim, we proceed as follows: we generate a (numerical) database, consisting of a set of concentration and velocity fields. Then, we use these concentration fields as input for the CVR algorithm. Finally, the velocity fields reconstructed with CVR are compared with those obtained numerically. We will first describe the numerical model used for the simulations and then compare the results obtained in terms of velocity distributions.

We assume that the flow field $\mathbf {u}$ is incompressible and described by the Darcy law, i.e., it obeys the continuity (3) and the Darcy Eq. (5). In the absence of dispersion and assuming a constant diffusion coefficient, solute concentration is controlled by the advection–diffusion equation (Nield et al. 2013, p. 42):

$$\begin{aligned} \frac{\partial C}{\partial t}+\mathbf {u}\cdot \nabla C=D\nabla ^{2} C{ .} \end{aligned}$$

(12)

The idea behind this model is that the time scale of molecular diffusion is much smaller than the convective time scale or, in other terms, the effect of inertia in negligible. This assumption is reasonable particularly in the instance of flows in narrow gaps. We consider a rectangular domain (height $H=50\,\hbox {mm}$ and length $L=100\,\hbox {mm})$. Since we want to validate the CVR algorithm in the instance of Darcy flow $(\epsilon ^{2}{\textit{Ra}}\rightarrow 0)$, we assumed a very narrow gap, $b=0.03$ mm, to estimate the permeability value used in the two-dimensional simulation. We discretise Eqs. (5) and (12) numerically in space using second-order finite volume and sixth-order compact finite-difference schemes, respectively, in a uniform domain consisting of $400\times 200$ nodes. Discretisation in time is achieved using an explicit third-order Runge–Kutta scheme. The code has been adapted from Hidalgo et al. (2012, 2015), where further numerical details are available. The boundary conditions are the same shown in Fig. 5, and are here summarised. For the solute transport equation (12), concentration is fixed along the top boundary ($C=C_{s}$) and no-flux condition is assumed at the side and lower boundaries ($\partial _{\mathbf {n} }C=0$). For the flow field, no penetration is applied along all the boundaries ($\mathbf {u}\cdot \mathbf {n}=0$). The domain is initially filled with pure water, i.e., $C(x,z<H,t=0)=0$, and characterised by a velocity $\mathbf {u}=(u,w)=0$. The concentration field is initially perturbed with a random noise of amplitude $10^{-3}C_s$. We remark here that the initial perturbation influences only the onset of convection, whereas the late convective dynamics is independent of the perturbation amplitude (Slim 2014; De Paoli et al. 2017).

Fig. 6

Validation of CVR in Darcy flow regime. a Contour map of the concentration distribution obtained from the numerical simulation. The iso-contours of concentration (black lines) are also shown. A small portion of the domain, indicated by the black rectangle, is considered in panels (b)–(c). Detail of the velocity field (vectors) obtained from numerical simulations (b) and reconstructed via CVR algorithm (c). The scale for the velocity vectors is also reported. The horizontal and vertical velocities measured along the two horizontal lines in panel (a) are reported in panels (d)–(g): Panels (d)–(e) and panels (f)–(g) correspond to $z=45~\hbox {mm}$ [solid line in panel (a)] and $z=25~\hbox {mm}$ [dashed line in panel (a)] respectively. The CVR results (symbols) are in excellent agreement with the numerical data (solid lines)

A snapshot of the concentration field obtained numerically is reported in Fig. 6a, where the concentration distribution is shown together with the iso-contours of concentration. Numerical data are used for the validation as follows. The horizontal concentration gradient in Eq. (6) is computed from the numerical concentration field C(x, z). Following the algorithm described in Sect. 3.1, we obtain the stream function and reconstruct the velocity field. A comparison of the velocity fields obtained with numerical simulations and CVR is provided in Fig. 6b–c, where the fingers, identified by the black rectangle in Fig. 6a, are considered. We observe qualitatively that the flow field (black vectors) is similar in both Fig. 6b, c, and the flow structures are well captured by the CVR. The core of the finger is characterised by a strong downward velocity and, due to continuity, the flow in the region between two fingers has a remarkable upward velocity that supplies the upper layer with fresh fluid. A more quantitative validation is proposed in Fig. 6d–g: the horizontal (u) and vertical (w) velocity components are measured along the solid ($z=45~\hbox {mm}$) and the dashed lines in Fig. 6a ($z=25~\hbox {mm}$), respectively, shown in Fig. 6d–e and 6f–g. We observe that, in correspondence of both values of z considered, the results obtained via CVR are in excellent agreement with the numerical data. The reconstructed velocity is compared with the numerical results for a long time range, and further results can be found in the electronic supplementary material available online (Movie 1).

This validation refers to the case of Darcy flows $\left( \epsilon ^{2}{\textit{Ra}}\rightarrow 0\right)$, which represents an ideal situation that can be attained in the instance of infinitesimal gap thickness, b, and low-density differences, $\Delta \rho$. However, the dynamics in the case of Hele–Shaw flows is controlled by the effect of the gap-induced solute dispersion, and some corrections to the Darcy equation should be taken into account (see Appendix A and electronic supplementary material).

3.3 Quantification of uncertainties

We analyse here the uncertainty on the values of the velocity reconstructed. The uncertainty values are obtained from propagation of error analysis (Taylor 1997). For the velocity components reconstructed as in Eq. (5), the uncertainty is given by:

$$\begin{aligned} \frac{\delta u}{u}= & {} \biggl [\left( 2\frac{\delta b}{b} \right) ^{2}+\left( \frac{\delta C}{C}\right) ^{2}+\left( \frac{\delta \mu }{\mu } \right) ^{2}\nonumber \\&+\left( \frac{\delta \Delta \rho }{\Delta \rho } \right) ^{2}+\left( \frac{\delta C_{s}}{C_{s}} \right) ^{2}+\left( \frac{\delta g}{g} \right) ^{2}\biggr ]^{1/2}\text { ,} \end{aligned}$$

(13)

where $\delta x_{i}$ indicates the uncertainty on the quantity $x_{i}$. We assume that the properties $\mu$, $\Delta \rho$, $C_{s}$, and g are estimated with constant relative uncertainty equal to $1\%$. We consider the accuracy on the gap thickness $\delta b=10\mu$m. All the uncertainties introduced are constant and independent of the value of concentration measured, but the uncertainty on the concentration itself, $\delta C$, varies with the value of local concentration. In particular, since $C=\rho \omega$, we estimate the uncertainty on C as:

$$\begin{aligned} \frac{\delta C}{C}=\sqrt{\left( \frac{\delta \rho }{\rho } \right) ^{2}+\left( \frac{\delta \omega }{\omega } \right) ^{2}}\text { .} \end{aligned}$$

(14)

The mass fraction is obtained from the transmitted light intensity ratio, $I(\omega )/I(0)$, as reported in Fig. 2b. We compute $\delta [I(\omega )/I(0)]=\sigma$, where $\sigma$ indicates the standard deviation of $I(\omega )/I(0)$, measured over ten images in correspondence of each mass fraction. Then, the uncertainty $\delta \omega$ is estimated from $\delta [I(\omega )/I(0)]$ (Fig. 2b). Although the uncertainty $\delta \omega$ is very limited for high-light intensities, we observed that the corresponding relative uncertainty is independent of $I(\omega )/I(0)$, and it is $\delta \omega /\omega \approx 0.08$. We approximate the uncertainty on the density measurement as $\delta \rho /\rho \approx \delta \rho (0)/\rho (0)\le 0.02\%$. Therefore, the contribution of the uncertainty of the mass fraction to the uncertainty of the concentration measurement is dominant, and we assume $\delta C/ C\approx \delta \omega /\omega$. As a result, due to Eq. (13), the uncertainty on the velocity measurement for a cell of thickness $b=0.30~\hbox {mm}$ is $\delta u / u =0.1$.

4 Example application: tip splitting dynamics

The mechanism by which the tip of a finger separates in two branches is defined as tip splitting (Malhotra et al. 2015). This phenomenon is currently subject of interest for the implications which it can bear in porous media flows (Amooie et al. 2018), e.g., the enhancement of the dissolution efficiency in buoyancy-driven flows (Jha et al. 2011). In the context of miscible fluids, the first detailed study on this phenomenon dates back to the seminal work of Tan and Homsy (1988). With the aid of numerical simulations, they have shown that once a finger is sufficiently large, the concentration gradient across the finger front becomes steep, as a result of the stretching produced by the flow. In turn, the tip of the finger becomes unstable and splits. When the finger splits into two parts and both new fingers continue to grow, the mechanism is defined as even tip splitting. When after the splitting, one finger shields the growth of the other one, the mechanism is defined as uneven tip splitting (Malhotra et al. 2015). More complicated scenarios occur when the branches are more than two. The tip splitting phenomenon is the result of complex non-linear interactions between the fingers, which makes predictions on its evolution hard to obtain. As a result, a variety of different flow configurations has been studied in the last decades. Numerical works (De Wit and Homsy 1997; De Wit 2004) report that variations in the permeability of the medium and domain size have an influence on the splitting process of the fingers. The presence of reactive components in the mixture will also foster the tip splitting process, via a mechanisms called chemically induced tip splitting (De Wit and Homsy 1999). A similar behaviour is observed when the permeability distribution of the medium is anisotropic (Kawaguchi et al. 2001). The situation becomes even more complicated when three-dimensional flow dynamics are taken into account.

In the numerical work of Oliveira and Meiburg (2011), the full three-dimensional flow in a Hele–Shaw cell is investigated in detail for the first time. They observed that fingers may split also longitudinally, and not only at their tip, giving rise to a different phenomenon defined as inner splitting. This result is also supported by the previous experimental work of Wooding (1969). Recently, the first three-dimensional experimental observation of tip splitting in porous media has been reported by Suekane et al. (2017). The authors preformed concentration measurements in a cylindrical porous medium, and observed that the splitting mechanism is similar to that observed in the two-dimensional case. However, an accurate experimental characterisation of the flow around the fingers during the tip splitting is still missing and we aim precisely at this gap. In the following, we will first investigate the tip splitting phenomenon via detailed analyses of the concentration and vorticity fields (Sect. 4.1). Then, using this information, we briefly present two different splitting dynamics observed in our experiments (Sects. 4.2–4.3).

Fig. 7

Tip splitting of a finger (gap width $b=0.30~\hbox {mm}$). Time-dependent evolution of solute concentration distribution (top panels) and vorticity contours (bottom panels) are here reported. Time is indicated in each panel. The vorticity contours are coloured according to the value of the vorticity: positive (negative) values of vorticity correspond to blue (red) contours and are associated with vortices rotating in clockwise (counter-clockwise) direction

4.1 Vortex formation and tip splitting

We consider an experiment characterised by $b=0.30~\hbox {mm}$. The concentration field measured experimentally is used to reconstruct the velocity field via CVR. In the present configuration, since $b=0.30~\hbox {mm}$, non-Darcy corrections have been applied (see Appendix A). The corresponding vorticity field, defined as $\varvec{{\Omega }}(x,z)=\nabla \times \mathbf {u}$, is then obtained and used for the analysis of the flow dynamics in correspondence of the finger tip. The evolution of the concentration field, C, and out-of-plane vorticity component, $\Omega _{y}$, are shown in Fig. 7, in top and bottom panels, respectively. The vorticity contours are coloured according to the value of the vorticity: positive (negative) values of vorticity correspond to blue (red) contours and are associated with vortices rotating in clockwise (counter-clockwise) direction. We observe that the evolution of the fingers is characterised by a complex dynamics. Initially, the finger spreads and the width of the structure increases (Fig. 7a, f). This process is driven by the flow field: The fluid is redirected away from the axis of the finger. As a result, the concentration gradient at the finger tip becomes steeper, preparing the ground for the instability to take place. Two pairs of counter-rotating vortices have the effect of stretching the structure. Moreover, small deviations from the vertical direction produce a symmetry breaking in the finger structure, making the finger unstable and promoting, eventually, the bifurcation of the tip (Fig. 7b, g). The two counter rotating vortices coexist for about 10 s (Fig. 7c, h to d, i), during which they grow in a nearly symmetric fashion. Finally, approximately 30 s after the tip splitting has started (Fig. 7e, j), the right-hand branch detaches from the main body due to a slightly unbalanced growth of the two fingers tips.

We conclude that the tip splitting observed here is controlled by the local velocity field. This observation is the first accurate simultaneous measurement of velocity and concentration field and is in excellent agreement with the previous numerical observations (Chen and Meiburg 1998).

4.2 Uneven tip splitting

We consider now an experiment again characterised by a cell gap $b=0.30~\hbox {mm}$, but we focused on a different phenomenon consisting of the differential growth of the finger branches. An example is shown in Fig. 8. The finger is initially symmetric and grows downward [8a,e]. However, when the splitting takes place [8b,f], a pair of counter rotating vortices forms at the finger tip. This time, their growth is strongly unbalanced from the beginning, and within few seconds, the left-hand finger is considerably larger than the right-hand counterpart [8c,g]. The amount of solute present in the two branches is also unbalanced: The larger finger mixes within surrounding fluid more slowly compared to the right-hand side branch. In this case, the combined action of flow and convective dissolution leads to the formation of a dominant finger that prevents the growth of the other branch. This phenomenon is also known as shielding (Tan and Homsy 1988).

Fig. 8

Uneven tip splitting of a finger (gap width $b=0.30~\hbox {mm}$). Time-dependent evolution of solute concentration distribution (top panels) and vorticity contours (bottom panels) are here reported. Time is indicated in each panel. The vorticity contours are coloured according to the value of the vorticity: positive (negative) values of vorticity correspond to blue (red) contours and are associated with vortices rotating in clockwise (counter-clockwise) direction

4.3 Post-merging tip splitting

Finally, we analyse the most recurrent phenomenon which we observed in our experiments, consisting of finger merging and subsequent separation. We consider the system shown in Fig. 9. Initially, two fingers merge and form the large finger reported in Fig. 9a, e. The merging process took place only at the boundaries of the fingers: the footprint of the core of the two initial fingers is still apparent and represented by the two high-concentration regions [dark distinct parts in Fig. 9b, f]. This heterogeneity in the solute distribution triggers the subsequent splitting (Fig. 9c, g) which, in this case, is nearly symmetric (even splitting). However, we observed that in our experiments, the uneven splitting is more likely to happen. It is also interesting to note the presence of a stagnation region in correspondence of the finger splitting point (located approximately at $(x,z)=(30,125)~\hbox {mm}$).

Fig. 9

Post-merging tip splitting (gap width $b=0.30~\hbox {mm}$). Time-dependent evolution of solute concentration distribution (top panels) and vorticity contours (bottom panels) are here reported. Time is indicated in each panel. The vorticity contours are coloured according to the value of the vorticity: positive (negative) values of vorticity correspond to blue (red) contours and are associated with vortices rotating in clockwise (counter-clockwise) direction

Fig. 10

Process of Concentration-based Velocity Reconstruction. a Small portion of an experimental high-resolution image ($b=0.30~\hbox {mm}$). The fluid is characterised by a purple colour, typical of the solute used $(\hbox {KMnO}_{4})$. b Contour map of the concentration field (obtained from light intensities) and streamlines (reconstructed via CVR). c–d Distributions of the horizontal and vertical velocity components. Time-dependent reconstruction is available in the electronic supplementary material (Movie 2)

5 Conclusions

We proposed a new method for the simultaneous measurement of concentration and velocity fields of Newtonian fluids in convective Hele-Shaw flows. We developed an algorithm, defined here Concentration-based Velocity Reconstruction (CVR), to obtain the gap-averaged velocity distribution from one single snapshot of the concentration field. An example is shown in Fig. 10, where the input data, consisting of a high-resolution experimental image (Fig. 10a), are used to reconstruct the concentration field (Fig. 10b). The flow stream function is computed (Fig. 10b) via CVR, and the horizontal and vertical gap-averaged velocity components are reconstructed (Fig. 10c–d). The system considered is controlled by the complex interplay of convection and diffusion. If solute diffusion across the gap dominates over vertical advection (e.g., when the cell gap is infinitesimal), the flow can be considered as two-dimensional and described by a Darcy model. When the cell gap is increased, the velocity gradient in wall-normal direction promotes solute dispersion (Taylor dispersion), making the flow to depart from the ideal two-dimensional Darcy behaviour. However, when the strength of vertical convection is equivalent to that of transverse diffusion, the gap-averaged flow can still be considered two-dimensional and it is defined as a Hele–Shaw flow (Letelier et al. 2019; De Paoli et al. 2020). The CVR algorithm has been validated with two-dimensional numerical simulations and with particle tracking velocimetry measurements (further information is available in the electronic supplementary material), in the instance of Darcy and the Hele–Shaw flows, respectively. When the gap width is large enough to allow the formation of more than one finger in transverse direction, the flow is three-dimensional and the CVR algorithm cannot be applied.

We employed the CVR algorithm to analyse the phenomenon of tip splitting in convective Hele–Shaw flows. This well-known problem received renovated attention due to the beneficial effects which it can bear in mixing in porous flows. Via high-resolution reconstruction of the flow velocity around the fingers, we are able to analyse experimentally, for the first time with this level of detail, the flow configuration and dynamics during the splitting process.

In this work, we applied the CVR to a buoyancy-driven Hele–Shaw flow in one-sided configuration. However, this method can be easily adapted to different systems with minor modifications (e.g., Rayleigh–Taylor flows in confined Hele–Shaw geometries, De Paoli et al. 2019). A further advantage of the present method consists of the possibility of performing velocity measurements at high resolution and on the whole domain, including the finger cores, where solute concentration is high and other methods based on particle tracking may fail. Finally, the absence of laser sources makes current approach a cost-effective and safe alternative to classical PIV/PTV methods. The need of a single camera and the possibility of relying on one picture for a detailed reconstruction of the solute concentration and velocity fields, reducing costs and complexity of the experimental facilities, represent further advantages of the reconstruction method proposed.

Appendix A: Velocity reconstruction of Hele–Shaw flows

Convective Hele–Shaw flows characterised by large values of the cell gap, b, are not homogeneous across the cell gap. However, when $\epsilon ^{2}{\textit{Ra}}\ll 1$, the system can still be considered two-dimensional. Letelier et al. (2019) derived analytically a two-dimensional flow model of this thin three-dimensional system. As a result of their mathematical derivation, based on a perturbed form of the Navier–Stokes equation, the effect of gap-induced solute dispersion appears as a correction of the Darcy model. The effect of these perturbative corrections has to be taken into account within the process of velocity reconstruction. The model has been derived analytically from Letelier et al. (2019) in the frame of natural convection and constant temperature difference between the top and bottom boundaries. In the context of mass transfer and in non-dimensional form, the momentum equation reads:

$$\begin{aligned} \epsilon ^2\frac{{\textit{Ra}}}{{\textit{Sc}}}\left[ \frac{6}{5}\frac{\partial \mathbf {u}}{\partial t}+\frac{54}{35}\left( \mathbf {u}\cdot \nabla \right) \mathbf {u}\right]= & {} -\nabla p-\mathbf {u}+C\hat{z}+\frac{6}{5}\epsilon ^2\nabla ^2\mathbf {u}\nonumber \\&-\frac{2}{35}\epsilon ^2{\textit{Ra}}\left( \mathbf {u}\cdot \nabla C\right) \hat{z} \text { ,} \end{aligned}$$

(15)

where ${\textit{Sc}}=\mu /\rho D$ is the Schmidt number, standing for the ratio of momentum to mass diffusivity, and $\hat{z}$ is the unit vector aligned with the vertical axis, z. The influence of inertia has been investigated in this specific configuration by De Paoli et al. (2020), who observed that it has a minor role and, therefore, the corresponding terms in the momentum equation can be neglected. Moreover, we have here that the Schmidt number is large, ${\textit{Sc}}\sim O(10^{3})$. As a result of these assumptions, Eq. (15) reduces to:

$$\begin{aligned}&\underbrace{\mathbf {u}=-\frac{b^{2}}{12\mu }\left[ \nabla P+\frac{\Delta \rho }{C_{s}}C\mathbf {g}\right] }_{\text {Darcy model}}\nonumber \\&\quad + \underbrace{\frac{b^{2}}{10}\nabla ^2\mathbf {u}-\frac{1}{210}\frac{b^{4}g\Delta \rho }{12C_{s}D\mu }(\mathbf {u}\cdot \nabla C)\hat{z} }_{\text {non-Darcy corrections}}\text { .} \end{aligned}$$

(16)

By taking the curl of Eq. (16), the momentum equation can be written in terms of stream function as:

$$\begin{aligned} \nabla ^{4} \psi -\frac{10}{b^2}\nabla ^{2}\psi= & {} -\beta \frac{\partial C}{\partial x}+\gamma \left( -\frac{\partial ^2 \psi }{\partial x^2}\frac{\partial C}{\partial z} \right. \nonumber \\&\left. -\frac{\partial \psi }{\partial x}\frac{\partial ^2 C}{\partial x\partial z}+\frac{\partial ^2\psi }{\partial x \partial z}\frac{\partial C}{\partial x}+\frac{\partial \psi }{\partial z}\frac{\partial ^2 C}{\partial x^2}\right) , \end{aligned}$$

(17)

where $\beta =10\alpha /b^2$ and $\gamma =\alpha /(21D)$. Numerical discretisation adopted for the Darcy flow [Eq. (6)] is used to solve Eq. (17). We remark here that Eq. (17) consists of a Helmholtz problem and the numerical solution is computationally more demanding than the corresponding Darcy problem. Moreover, due to the high order of the derivatives in Eq. (17), the solution is more sensitive to the input data. The same boundary conditions considered for the Darcy model and discussed in Sect. 3.1 are applied here. In the instance of Hele–Shaw flows, the CVR has been validated against Particle Tracking Velocimetry (PTV) measurements, showing both quantitative and qualitative agreement. Further information is available in the electronic supplementary material.

Fig. 11

Correlation factor r between the Darcy and the Hele–Shaw model, averaged from the beginning to the shutdown phase, for 18 experimental data (De Paoli et al. 2020). When $\epsilon ^{2}{\textit{Ra}}\le 0.01$, the results given by the CVR with the Darcy and the Hele–Shaw model are in excellent agreement. For larger values of $\epsilon ^{2}{\textit{Ra}}$, the correlation factor diminishes, indicating that the prediction given by the Darcy model will be less and less accurate when the role of convection becomes dominant

We analysed the results of the CVR algorithm in the instance of Hele–Shaw flows. In particular, we compared the predictions of the Darcy model, Eq. (4), with the Hele–Shaw model, Eq. (16), in different experimental conditions. To this aim, we used the database of De Paoli et al. (2020), consisting of 18 experiments split over three different gap thicknesses ($b=0.15$, 0.30 and 0.50 mm). Each test has been repeated three times, and the results reported represent the average of these three experiments. The flow field is reconstructed with the CVR algorithm using the Darcy model and the Hele–Shaw model. To evaluate the difference in the predictions given by the Darcy and the Hele–Shaw model, following Kähler et al. (2016), we use the Pearson’s correlation coefficient, r. We estimate, in each sample corresponding to time t, the space-average magnitude of the local fluid velocity:

$$\begin{aligned} \bar{V}(t)=\frac{1}{N_{x}N_{z}}\sum _{ij} V(i,j,t)\text { ,} \end{aligned}$$

(18)

with $V(i,j,t)=\sqrt{[u(i,j)]^{2}+w[(i,j)]^{2}}$ the local fluid velocity magnitude, i, j the pixel indices, and $N_{x}$, $N_{z}$ the image dimensions. Then, we compute the time-averaged correlation coefficient r between the velocity magnitude reconstructed in the Darcy (D) and in the Hele–Shaw (HS) cases as:

$$\begin{aligned} r=\frac{1}{T}\int _{0}^{T} \frac{\sum _{ij} \left( V_{D}-\bar{V}_{D}\right) \left( V_{HS}-\bar{V}_{HS}\right) }{ \sqrt{\sum _{ij} \left( V_{D}-\bar{V}_{D}\right) ^{2}}\sqrt{\sum _{ij} \left( V_{HS}-\bar{V}_{HS}\right) ^{2}} } \text { d}t\text { ,} \end{aligned}$$

(19)

where T is the time at which the domain saturation occurs (see electronic supplementary material for further details). Results are shown in Fig. 11 in terms of $r(\epsilon ^{2}{\textit{Ra}})$. When $\epsilon ^{2}{\textit{Ra}}\le 0.01$, the value of the correlation coefficient is close to 1, suggesting that the velocities obtained with the two models are strongly correlated and, therefore, we are in the presence of a Darcy flow. For larger values of $\epsilon ^{2}{\textit{Ra}}$, the value of the correlation coefficient r is lower, but still positive. This means that the velocities, $V_{D}$ and $V_{HS}$, lie on the same side of their respective means, $\bar{V}_{D}$ and $\bar{V}_{HS}$, but are less correlated compared to the cases when $\epsilon ^{2}{\textit{Ra}}\le 0.01$. The correlation factor diminishes when $\epsilon ^{2}{\textit{Ra}}$ is further increased: The prediction given by the Darcy model is less and less accurate and vertical advection has a more dominant role.