How non-Darcy effects influence scaling laws in Hele-Shaw convection experiments

We examine experimentally the influence of non-Darcy effects on convective dissolution in Hele-Shaw cells. We focus on buoyancy-driven convection, where the flow is controlled by the Rayleigh–Darcy number, $Ra$ , which measures the strength of convection compared to diffusion. The Hele-Shaw cell is suitable to mimic Darcy flows only under certain geometrical constraints, and a recent theoretical work (Letelier et al., J. Fluid Mech., vol. 864, 2019, pp. 746–767) demonstrated that a precise limit exists for the parameter $\unicode[STIX]{x1D716}^{2}Ra$ – $\unicode[STIX]{x1D716}\sim$ thickness-to-height ratio – beyond which the flow exhibits non-Darcy effects. In this work, we run experiments for solute convection in Rayleigh–Bénard-like configuration. We examine a wide range of the parameters space $(Ra,\unicode[STIX]{x1D716})$ and we clearly identify the application limits of Darcy flow assumptions. Besides confirming previous theoretical predictions, current results are of relevance in the context of porous media flows – which are often studied experimentally with Hele-Shaw set-ups. Using our original datasets, we have been able to explain and reconcile the discrepancies observed between scaling laws previously proposed for Rayleigh–Bénard-like experiments and simulations in similar contexts. Specifically, we attribute an important role to the parameter $\unicode[STIX]{x1D716}^{2}Ra$ , which clearly establishes thresholds beyond which Hele-Shaw experiment results are influenced by three-dimensional effects.


Introduction
Hele-Shaw cells are made by two parallel, transparent plates, which are separated by a gap that -ideally -is infinitesimal. Under this condition, any flow established inside the gap can reproduce a Stokes flow. This property makes Hele-Shaw cells widely used to visualise and analyse experimentally different types of flow instances. We focus here on the use of Hele-Shaw cells to investigate flows driven by the differential gravity force acting on fluids of different density, which are initially set in an unstable configuration. This type of experiment is particularly important for a number of applications, since it is simple and can, under certain assumptions, reproduce important features of Darcy flows in porous media. Specifically, we refer to the instance of fluid saturated porous media, in which a heavier fluid sits on top of a lighter fluid. In this instance (for certain density differences), the fluid flow is dominated by convection; this type of flow was examined by Saffman & Taylor (1958), who also introduced the concept that the flow inside a porous media can be experimentally reproduced employing a Hele-Shaw cell.
In this paper, we examine the process of solute convection in Rayleigh-Bénard-like flows in Hele-Shaw cells, with the following set-up: the cell is filled with pure fluid, i.e. solute is initially not present, and the concentration of solute is kept constant along the top boundary of the domain. Since the density of the fluid is larger at the top wall, the system is unstable and, with time, tiny instabilities of the density interface grow into finger-like structures, which eventually develop convective currents. This configuration is of crucial importance to investigate geophysical subsurface flows such as water contamination (LeBlanc 1984;Van Der Molen & Ommen 1988), petroleum migration (Simmons, Fenstemaker & Sharp Jr 2001), carbon sequestration (Huppert & Neufeld 2014;Emami-Meybodi et al. 2015) and sea ice formation (Wettlaufer, Worster & Huppert 1997;Feltham et al. 2006). The main dimensionless governing parameter of the flow is the Rayleigh-Darcy number (Ra), which stands for an inverse diffusivity, or for the relative strength of convection to diffusion (Slim 2014), and is defined as where, ρ * s is the density difference between the solute-saturated fluid and the pure fluid, g is the acceleration due to gravity, b * is the cell gap thickness, H * is the cell height, D is the molecular diffusion coefficient and µ is the dynamic viscosity. This system is analysed in terms of solute dissolution rate F(t), i.e. the amount of solute dissolved per unit of area and time. Based on the time evolution of F(t), three different flow phases can be defined: (i) the flow is initially controlled by diffusion; then (ii) an unstable layer builds up at the upper boundary and convective finger-like structures form and control the evolution of the system; and finally (iii) shutdown of convection takes place (Hewitt, Neufeld & Lister 2013;Slim 2014).
Two-dimensional Darcy simulations have shown that, during the convection dominated phase, the solute dissolution rate is independent of the Rayleigh-Darcy number (Pau et al. 2010;Hewitt et al. 2013;Slim 2014). Noteworthy though, experiments in Hele-Shaw cells have shown that the quantity F(t) exhibits weak dependence on Ra (Backhaus, Turitsyn & Ecke 2011;Tsai, Riesing & Stone 2013). Therefore, in an attempt to unravel this discrepancy, Hidalgo et al. (2012) investigated numerically the time-dependent dissolution process in fluids characterised by concentration-dependent viscosity and non-monotonic density-concentration curves; they concluded that the Ra-dependent character of the dissolution flux observed experimentally is not motivated by these two causes. Subsequent studies investigated the influence of the geometry of the domain on the features of the flow. In particular, with the aid of perturbation techniques and numerical simulations, Letelier, Mujica & Ortega (2019) investigated the effect of the combined action of the Rayleigh-Darcy number and the anisotropy ratio = b * / √ 12H * . In accordance with the value of the parameter 2 Ra, they have been able to identify three flow configurations: (i) Darcy regime ( 2 Ra → 0), where the flow is two-dimensional and well described by Darcy simulations; (ii) Hele-Shaw regime ( 2 Ra 1), where the flow is still two-dimensional but influenced by gap-induced dispersion; and (iii) three-dimensional regime ( 2 Ra 1), when the effects of the third dimension become non-negligible.  The object of this work is to provide experimental ground to these recent findings, with a further ambitious attempt of reconciling the different scaling laws of the solute dissolution rate available from the literature. To explore a wide range of the parameter space 2 Ra, we will use a Hele-Shaw cell with variable gap and variable height in a Rayleigh-Bénard-like set-up.

Methodology
We consider a rectangular domain initially filled with pure fluid. Solute concentration is kept constant at the top wall, whereas other boundaries are impermeable with respect to fluid and solute fluxes. This set-up is defined as Rayleigh-Bénard-like or a one-sided configuration (Hewitt et al. 2013;De Paoli, Zonta & Soldati 2016). For different Rayleigh-Darcy numbers, we measured the mean solute dissolution rate by keeping the same fluid and solute ( ρ * s , D and µ are constant), but varying the geometry of the cell (gap width, b * , and domain height, H * ).
The apparatus is sketched in figure 1(a) and consists of two parallel acrylic plates (thickness 30 mm), having width of 200 mm and height of 370 mm, separated by a gap b * ∈ [0.15; 1.00] mm. Except the top wall, fluid layer boundaries are defined by rubber seals, which confine the fluid between the plates and avoid leakages. The cell width is kept constant and equal to L * = 160 mm, whereas the cell height is varied using different gaskets so that H * ∈ [104; 343] mm. The material used for the gaskets is a high-quality impermeable rubber (Klinger-Sil C-4400) worked with highprecision CNC machines. The fluid domain, with explicit indication of the dimensions, coordinate system and boundary conditions, is shown in figure 1(b). The plates are held in place with the aid of screws. To obtain a desired value of gap thickness and to ensure its uniformity over the cell, metal shims are placed in the gap and the same torque is applied to all the screws. Along the top wall lies a steel mesh (40 µm grid size) which contains the dye powder.
Jafari- Raad & Hassanzadeh (2015) observed that some discrepancies between experiments and simulations may occur due to the nature of the fluids adopted. For instance, Methanol and Ethylene-Glycol (known as MEG) and Propylene-Glycol (PPG), often used in experimental studies, are characterised by a non-monotonic variation of the fluid density with respect to solute concentration. In contrast, most of the numerical studies considered a much simpler linear dependency of density and concentration. To neglect the effect of the fluid properties, we chose two fluids marked by a linear density-concentration profile: an aqueous solution of potassium permanganate (KMnO 4 ) to mimic the heavy fluid (high solute concentration) and water to mimic the light fluid (low solute concentration). Viscosity and diffusion are assumed constant and equal to µ = 9.2 × 10 −4 Pa s and D = 1.65 × 10 −9 m 2 s −1 , respectively (Slim et al. 2013). With these fluids, a linear behaviour of the mixture density ρ * as a function of the solute concentration C * is obtained, i.e.
where, ρ * s is the density of the saturated solution, C * s = 48 kg m −3 is the effective saturation value of concentration and ρ * s is the density difference between a saturated solution and pure water. A discussion on the mixture density as a function of the solute concentration can be found in appendix A. Since water density is sensitive to temperature, although we run experiments in the temperature range 20-25 • C, for the sake of precision ρ * s and ρ * s are estimated in each experiment by the correlations of Novotný & Söhnel (1988). With this set-up, we obtained Ra ∈ [4.8 × 10 4 ; 7.0 × 10 6 ]. We reconstruct solute concentration fields from light intensity maps (Slim et al. 2013;Ching, Chen & Tsai 2017). The accuracy of the concentration reconstruction process is sensitive to the local value of the mass fraction (Slim et al. 2013); in particular, the error may be significant (a few per cent) in regions of the domain that are characterised by high values of solute concentration. However, regions of high solute concentration are limited to the top of the domain, with an overall extension of less than 5 %, similar to previous experimental campaigns (Slim et al. 2013;Alipour & De Paoli 2019). Grains of potassium permanganate (diameter greater than 200 µm) are poured onto the grid. Afterwards, water is injected from two channels located at the bottom of the cell with the aid of a syringe pump. The fluid level is increased up to the upper boundary and the pump is shut down. This moment is considered as the beginning of the experiment (t * = 0). After water injection, the solute dissolves and a fluid layer denser than water forms below the grid: the heavy mixture layer thickens, becomes unstable and the finger-like structures form (Slim 2014). We use a Canon 1300D camera equipped with 18-55 mm lenses to record the evolution of the flow, with resolution and acquisition rate corresponding to 3456 × 5184 pixel and 1 f.p.s., respectively. The system is illuminated from the back side with a tuneable LED-based system: the tension applied to the LEDs is adjusted to maximise the sensitivity of the apparatus. Since the gap thickness has a significant effect on the colour of the mixture, the calibration is repeated for each value of b * . Given the light intensity distribution (camera images), the concentration field is finally reconstructed (Slim et al. 2013). Experimental results are the average of three experiments in the same configuration. A summary of all experimental parameters investigated can be found in table 1.
To present the results in dimensionless form, we set the buoyancy velocity as the reference velocity scale, whereas lengths and time are scaled with L * = D/W * and L * /W * , respectively. Concentration is made dimensionless with respect to C * s (for further details on the dimensionless set of variables, see Slim (2014)   of the superscript * is used here to refer to dimensionless variables. The most relevant global observable for the present system is the averaged dimensionless dissolution rate F(t). This observable, which represents the amount of solute dissolved per unit of area and time, is customarily computed in numerical simulations as F(t) = 1/L L 0 ∂ z C(x, z)| z=H dx. However, in experiments, the concentration gradients, especially at the top of the domain, are very sensitive to the concentration values as reconstructed from the pixels' intensity. Therefore, we compute F(t) starting from the cumulative indicator represented by the dimensionless mass of solute contained in the domain per unit of depth, m(t) = L 0 H 0 C dx dz, obtaining (Ching et al. 2017) Numerical and experimental investigations have shown that, with respect to the dissolution process, three regimes may be identified: diffusion-driven, convectiondominated (or constant flux) and shutdown of convection. We refer to Slim (2014) for a thorough description of the entire dissolution dynamics.

Results and discussion
In this section we will assess potential consequences of the assumption of Darcy flow by examining the role of the cell width when it is used as the parameter to increase the Rayleigh-Darcy number. However, it is instrumental in our analysis to examine the behaviour of dissolution flux in time for cases in which the Darcy flow assumptions can be safely applied (we will demonstrate this later). For experimental convenience, we will examine the solute flux behaviour for different Rayleigh-Darcy numbers: for this assessment, we increase Ra by increasing the domain height H * , as from (1.1), but fixing the gap thickness to b * = 0.30 mm. The evolution of the dimensionless time-dependent dissolution rate, F(t), computed as in (2.3), is shown in figure 2. The line at Ra = 6.3 × 10 5 is black as a reference for later discussion. It is apparent from the plateaux exhibited by the different experiments that a constant flux regime is found and it is later followed by a shutdown regime. The value of the plateau is similar for all the different cases and appears independent of Ra. sharp concentration gradients (highlighted by the pronounced colour gradients) are noticeable, we observe from figure 3(b-d) that increasing Ra via increasing the gap width corresponds to a smoothing of the concentration gradients, and fingers appear progressively less and less coherent. Responsible for this smoothing of the concentration gradients is the shear (or hydrodynamic) dispersion (Taylor 1953), which occurs when a velocity profile advects scalar species at different streamwise velocity. Increasing the gap width corresponds to a decrease of the resistance experienced by the fluid flowing down the cell and, for the same driving force (i.e. density difference) we should expect a larger solute flux, since the convective velocity W * scales with (b * ) 2 , see (2.2). However, since the hydrodynamic dispersion varies with the velocity gradient of the Poiseuille flow profile, which scales as . Larger gap thicknesses correspond to stronger effects of dispersion. When gap thickness is increased, the shape of the fingers is less defined and the concentration gradients across the interface of the fingers reduce.
increasing the gap width produces an interplay between convective velocity and hydrodynamic dispersion; in this way, the flow in the cell can be influenced by the three-dimensional spurious effects, which can prevent the application of the Darcy flow hypothesis to analyse the flow features. This will be further discussed later in this section.
To fully appreciate the important role of hydrodynamic dispersion on this type of experiment, we present in figure 4 a global view of all our measurements of F , which is the time average of the solute flux F(t) during the constant flux regime. We underline here that Letelier et al. (2019) used the scalar dissipation rate rather than F , but as discussed by De Paoli et al. (2019), for the present configuration, F is equivalent to the mean scalar dissipation rate. In figure 4(a), we present F as a function of Ra, which is varied both by changing H * and b * , see (1.1). The figure shows five datasets, corresponding to five values of b * . We can observe that data corresponding to small values of b * ( 0.50 mm) have a distinctly different behaviour from those obtained for the two larger values of b * . In particular, for b * = 0.80 mm and 1.00 mm, corresponding to a decrease of resistance to the flow, there is an important decrease of F , while hypothetical two-dimensional simulations would give a constant F . The case of the reported strong reduction of F in our experiments may be interpreted as the combined action of spurious transversal solute fluxes in the Hele-Shaw cell, which produce the dispersion effects as discussed by Letelier et al. (2019), and the transition towards to the three-dimensional flow regime. Numerical simulations can replicate a perfect two-dimensional situation and therefore will be used here to benchmark and analyse our data in the limit of small thicknesses. Numerical simulations are commonly analysed in terms of the Sherwood number, Sh = F Ra which, in the instance of heat transfer problems, is called the Nusselt number. In the papers by Slim (2014) and Wen et al. (2018a), a linear scaling was found for the Sherwood number as Sh = βRa α with α = 1 and β = 0.017. Similarly, a linear scaling is obtained by Hidalgo et al. (2012) and Amooie et al. (2018). With the aim of emphasising the effect of the experimental parameters on Sh, in figure 4(b) we show the behaviour of the Sherwood number as a function of Ra in the following way: we grouped the data into subdatasets of six experiments obtained for the same gap, and for each subdataset, Ra was varied by changing only the domain height H * . In figure 4(b), the scaling exponents found via data fitting and limited to each subdataset are also reported: the value of the exponents is in very good agreement with those found in the literature, which are summarised in table 2. To find best fit power laws, a regression model is applied on the data in the form log Sh = α log Ra + log β.
If we compare now our situation to previous experimental campaigns, we observe that Backhaus et al. (2011) and Tsai et al. (2013) report α = 0.80 and 0.76, respectively. In both papers α was identified as a best fit exponent on the entire dataset, but while Tsai et al. (2013) varied b * only, Backhaus et al. (2011) varied both b * and H * . To emphasise clearly the effect of changing b * and H * on our data, in figure 5 we show the behaviour F as a function of Ra and we find the best fit exponents α on different subdatasets which, this time, group data points with the same H * . We plot F rather than Sh because effects emerge in a more evident manner, and therefore we now look for the scaling   2. Summary of scaling exponents, α of equation Sh = βRa α , as found in previous numerical and experimental works. We consider here only experiments performed in Hele-Shaw cells and two-dimensional Darcy simulations in homogeneous and isotropic porous media. In the experimental works, an aqueous solution of PPG has been used. In the numerical works, different discretisation techniques have been adopted. The range of Rayleigh-Darcy numbers considered is also reported.
of α found in previous experimental works includes also some effect of the spurious transversal currents.
To analyse in detail, and also in an effort to quantify the influence of the geometry of the Hele-Shaw cell, which is crucial for the evolution of the flow, we use the cell anisotropy ratio, = b * / √ 12H * , which scales as the ratio of the characteristic time of transverse diffusion, (b * ) 2 /D, to the longitudinal advection time, H * /ŵ * , beinĝ w * the vertical velocity averaged across the gap b * (see Bae, Beta & Bodenschatz 2009;Kirby 2010). The behaviour of F as a function of the anisotropy ratio Results are grouped by regime as Hele-Shaw and three-dimensional regime. Here, (b) F is shown as a function of 2 Ra. In this case, the regime classification is more evident. A drop of the dissolution rate is observed in correspondence of 2 Ra = 1, where the transition from Hele-Shaw flow to three-dimensional flow occurs. The Darcy regime, which represents a theoretical limit for the experiments and corresponds to 2 Ra → 0, is also indicated.
is shown in figure 6(a). As discussed by Letelier et al. (2019), a Hele-Shaw cell behaves as a two-dimensional domain only in the limit of → 0, which is also the limit for Darcy flow. Indeed, an increase of the gap width produces an increase of vertical velocity, while the characteristic diffusion time scale remains unchanged. Since diffusion is responsible for the homogenisation of the solute distribution in the wall-normal direction, we cannot assume that the solute concentration is uniform across the gap. An increase of implies the occurrence of transverse effects which will lead the behaviour of the cell to Hele-Shaw regime and further to three-dimensional regime as proposed by Letelier et al. (2019). The application of their classification to our data emerges very clearly from figure 6(b), where we plot F as a function of the dimensionless number 2 Ra: if 2 Ra is small, the cell behaviour is closer to the Darcy flow behaviour, whereas if 2 Ra is large the flow exhibits three-dimensional effects which cannot be described with the Darcy equations, and can cause discrepancies of the scaling exponents. This is perhaps rather obvious, but it is here quantified precisely via experiments for the first time. We observe a drop of F after the threshold 2 Ra = 1, which represents the threshold value identified theoretically for the transition to the three-dimensional flow (Letelier et al. 2019). Our current dataset does not allow us to investigate whether there is a sharp or a smooth transition in correspondence of 2 Ra = 1. We can, however, confirm that the transition to the three-dimensional regime occurs for the subdatasets b * = 0.80 and 1.00 mm, corresponding to 2 Ra 1. These results provide a further assessment of the gap-induced dispersion effects, which represent one of the possible differences in the flow behaviour observed numerically and experimentally. Clearly, a transition of the cell behaviour from Hele-Shaw to three-dimensional regime, will require using the full Navier-Stokes equation set rather than a Darcy model (plus corrective terms) to analyse the flow. From the experimental viewpoint, a satisfactory analysis of full three-dimensional regime would require further observations in the out-of-plane direction. Phenomena like multiple finger in the thickness of the cell could not, otherwise, be accurately quantified.
The effect of inertia in the Hele-Shaw regime was examined theoretically by Letelier et al. (2019), who proposed to add extra terms to the Darcy model. In our experiments, these extra terms scale as 2 Ra/Sc, where Sc is the Schmidt number, i.e. the ratio of momentum to mass diffusivity, and is defined here as Sc = µ/ρD. We have that in our experiments Sc ∼ O(10 3 ). Since inertial terms are multiplied by the coefficient 2 Ra/Sc (Letelier et al. 2019), we can neglect these terms when 2 Ra 1, which in our experimental dataset corresponds to b * 0.50 mm. Our aim here is to evaluate the Reynolds number Re to confirm the predictions of Letelier et al. (2019). In this context, we assume as velocity scale w * f , defined as the mean value during the constant flux regime of the average fingertip velocity w * f , as sketched in figure 7(a), and we compute the Reynolds number Re as follows: To find w * f we use the horizontally averaged concentration profile We customarily set the location of the average fingertip as the vertical coordinate z * of the cell where C * (z * , t * )/C * s = 1/50, (3.4) as illustrated in figure 7(b). In figure 7(c) we report the values of Re for our dataset. We can observe that the Reynolds number Re increases for increasing the cell gap width, b * . The smallest gap-width experiments (b * = 0.15 mm) give a value of Re lower than 10 −1 , which is a good experimental approximation of Darcy flow conditions. For Re 1 (i.e. b * 0.50 mm) inertial effects are of the same order of viscous effects. However, we have shown via detailed analysis of the dissolution rate that these experiments are in the Hele-Shaw regime. This indicates that the Reynolds number is not sufficient for the evaluation of inertia, and a detailed analysis of the local flow velocities may be required. Finally, for the larger gap-width experiments Re 1 (i.e. b * 0.80 mm), where three-dimensional effects are important, we can observe that also inertial effects are becoming important, with Re value of the order of 10.

Conclusions
In this work, we investigate experimentally the problem of solute convection in Hele-Shaw cells in Rayleigh-Bénard-like configuration. The solute dynamics is governed by the complex interplay of convection, diffusion and dispersion. The flow is controlled by two dimensionless parameters: the Rayleigh-Darcy number, Ra, which measures the relative importance of convection and diffusion, and the anisotropy ratio, , proportional to the cell thickness-to-height ratio. On the basis of accurate measurements of the concentration field, we examined the occurrence of non-Darcy effects and their influence on the dissolution dynamics. Following the classification proposed by Letelier et al. (2019), based on the value of the dimensionless parameter 2 Ra, we can experimentally confirm the existence of three different flow regimes: (i) Darcy regime, for 2 Ra → 0, the flow is two-dimensional and is well described by a Darcy model; (ii) Hele-Shaw regime, the effect of non-Darcy terms is crucial for 0 < 2 Ra < 1, when dispersion effects are responsible for the Ra-dependent behaviour of the dissolution rate; (iii) three-dimensional regime, at last, for 2 Ra 1, the flow has a three-dimensional character. On the basis of this analysis, we have also been able to reconcile the Ra-dependent behaviour of the mean dissolution rate reportedly observed in previous experiments, and we have been able to attribute this dependency on Ra to spurious three-dimensional effects.
Present results can also have implications for porous media flows. To investigate the effect of a porous structure on solute convection, more complicated experiments use Hele-Shaw cells: the porous matrix is usually mimicked by glass beads, which fill the cell gap that in this case is of the order of few centimetres. Experimental investigations in Rayleigh-Bénard-like configuration indicate that, as well as in the Hele-Shaw apparatus, also in this case there is a regime dominated by convection in which the solute dissolution rate is constant, and the time-averaged dissolution rate scales as ∼Ra α−1 , with 0 < α < 1 (Neufeld et al. 2010). In a recent work, Liang et al. (2018) investigated experimentally the effect of solute redistribution induced by the beads. Tortuosity of the flow paths as well as friction with the solid surface of the porous matrix makes the fluid and the solute advected by the fluid follow sinuous pathways. This process, defined as mechanical dispersion, causes additional mechanical mixing and dilution effects and was identified as responsible for the Ra-dependent behaviour of the dissolution rate (Liang et al. 2018). This conclusion is also supported by simulations performed including the effect of mechanical dispersion (Wen et al. 2018b). The effect of mechanical dispersion in porous media shares some similarities with the Taylor hydrodynamic dispersion induced by the presence of the walls in Hele-Shaw cells. Although different in nature (e.g. strict two-dimensionality of the Hele-Shaw cell and tortuosity of flow paths in porous media), the Hele-Shaw cell is a good analogy for flows in homogenous and isotropic porous media provided that 2 Ra → 0, and can be used for the verification of closure models for Darcy flows (Nield & Bejan 2013).
where A 0 (T * ), A 1 (T * ) and A 2 (T * ) are functions of the temperature T * . In figure 8, the value of the density difference between solution (concentration C * ) and pure water is shown as a function of the solute concentration (symbols) for T * = 23 • C. The linear function corresponding to (2.1) (solid line) fits very well the correlation proposed by Novotný & Söhnel (1988) in the range of values of interest.