Role of solutal free convection on interdiffusion in a horizontal microfluidic channel
RRole of solutal free convection on interdiffusionin a horizontal microfluidic channel
Jean-Baptiste Salmon
CNRS, Solvay, LOF, UMR 5258, Univ. Bordeaux, F-33600 Pessac, France.
Laurent Soucasse
Laboratoire EM2C, CNRS, CentraleSup´elec,Universit´e Paris-Saclay, Grande Voie des Vignes,92295 Chatenay-Malabry cedex, France.
Fr´ed´eric Doumenc
Universit´e Paris-Saclay, CNRS, FAST, 91405, Orsay, France, andSorbonne Universit´e, UFR 919, 4 place Jussieu, F-75252, Paris Cedex 05, France. (Dated: February 19, 2021)
Abstract
We theoretically investigate the role of solutal free convection on the diffusion of a buoyant solute at themicrofluidic scales, (cid:39) µ m. We first consider a horizontal microfluidic slit, one half of which initiallyfilled with a binary solution (solute and solvent), and the other half with pure solvent. The buoyant forcesgenerate a gravity current that couples to the diffusion of the solute. We perform numerical resolutions ofthe 2D model describing the transport of the solute in the slit. This study allows us to highlight differentregimes as a function of a single parameter, the Rayleigh number Ra which compares gravity-inducedadvection to solute diffusion. We then derive asymptotic analytical solutions to quantify the width of themixing zone as a function of time in each regime and establish a diagram that makes it possible to identifythe range of Ra and times for which buoyancy does not impact diffusion. In a second step, we presentnumerical resolutions of the same model but for a 3D microfluidic channel with a square cross-section. Weobserve the same regimes as in the 2D case, and focus on the dispersion regime at long time scales. We thenderive the expression of the 1D dispersion coefficient for a channel with a rectangular section, and analysethe role of the transverse flow in the particular case of a square section. Finally, we show that the impactof this transverse flow on the solute transport can be neglected for most of the microfluidic experimentalconfigurations. a r X i v : . [ phy s i c s . f l u - dyn ] F e b . INTRODUCTION Microfluidics refers to a wide range of technological tools for manipulating liquids in microfab-ricated networks of channels with cross-sectional dimensions ranging from a few microns to a fewhundred microns. Applications of this technology are numerous and diverse, from high through-put miniaturized bioassays to fundamental studies in physical-chemistry, see Refs. [1–3] for somereviews. The small scales of microfluidic technologies allow to study numerous processes while Z X Y g (a) Z X W gravity current g HLH log T l og W b u o y a n c y e ff e c t s (b) U B W ∼ √ D T FIG. 1. (a) Schematic diagram of the diffusion experiment in a microfluidic channel of height H andwidth L . The red arrows represent the gravity current U B induced by the difference in density. The colorcode indicates the concentration of the solute. (b) Width of the mixing zone W vs. time T . The dashed lineis the case for which buoyancy does not affect the spreading of the solute. controlling finely all transport phenomena (mass, momentum, energy) [4]. In particular, the role ofbuoyancy has been mentioned by Squires and Quake’s in their review on the physics of fluids at thenanoliter scale and quantified using scaling arguments as explained below [5]. In this reference,the authors considered the situation illustrated in Fig. 1(a): a microfluidic channel of height H ,one half of which is initially filled with a binary solution, solvent and solute at concentration Φ i ,the other half only by the solvent. With no difference in density, solute and solvent interdiffuse,and the width of the mixing zone evolves as W ∼ √ DT where D is the diffusion coefficient of themixture and T the time. Now assuming that the density evolves linearly with the concentration insolute, i.e.: ρ = ρ ( + β Φ ) , (1)2here ρ is the density of the solvent, buoyancy induces a gravitational current, the solution flow-ing under the solvent for β >
0, see Fig. 1(a). The magnitude of the gravity current U B can beestimated from a balance between buoyant forces ∼ ρ β Φ i g and viscous forces ∼ ρ ν U B / H ,leading to the velocity scale: U B ∼ β Φ i gH ν , (2)where ν is the kinematic viscosity of the mixture [5]. Note that such a flow exists whatever theheight H of the channel as the density gradient is orthogonal to the gravity field g . The impact ofthis flow on the solute transport can be determined using the Rayleigh number:Ra = β Φ i gH ν D ∼ U B HD , (3)comparing advection to diffusion. Ra is similar to the P´eclet number describing mass transport inforced convection problems, but with the important difference that the solute in the present caseis not passive , since the solute itself is the source of the flow [5]. Furthermore, a balance betweenviscous forces ∼ ρ ν U B / H and inertial forces ∼ ρ U / H leads to the definition of the Grashofnumber, similarly to the Reynolds number for forced convection:Gr = β Φ i gH ν = RaSc , (4)where Sc = ν / D is the Schmidt number. As Sc ≥ for most liquid mixtures, one has thusGr (cid:28) Ra and viscous dissipation a priori dominates inertia within the gravity current in mostmicrofluidic applications [5].Many groups have reported gravity-driven currents in microfluidic experiments for which den-sity gradients are imposed either using membranes [6] or transverse mixing between coflowingmiscible liquids [7, 8]. In a different context, many groups also reported such flows when densitygradients are induced by the evaporation of a liquid mixture in a confined geometry ( H = µ m) such as sessile drops [9–11], confined drops [12–15], or micro-capillaries [16, 17].The impact of buoyancy on the solute transport is not always mentioned in such works, and mostgroups consider that solutal free convection plays little role at the microfluidic scales althoughexperimental configurations with density gradients are ubiquitous in applications.To illustrate this point, let us consider interdiffusion between an aqueous NaCl solution ata concentration of 1 M and pure water. For such a mixture, D (cid:39) . × − m /s, the dif-ference in density is (cid:39)
38 kg/m leading to β Φ i (cid:39) . × − , and the kinematic viscosity is3 (cid:39) − m /s [17]. For a microfluidic channel of height H = µ m, one finds Ra (cid:39) .
9, butdue to the scaling Ra ∝ H , the Rayleigh number increases to Ra (cid:39)
230 for H = µ m, andeven up to Ra (cid:39) . × for H = µ m. This numerical application illustrates the importanceof going beyond the scaling laws presented above to quantitatively predict the range of Rayleighnumbers for which free convection plays only a minor role in a microfluidic configuration. Further-more, although these buoyancy-induced flows could have little impact on the solute concentrationgradients that generate them, they still do exist, and are able to effectively disperse less mobilespecies in the fluid mixture, such as macromolecules or colloids [8]. These buoyancy-driven flowsmay also have an influence in protein crystallization experiments [18–21], for evaluating colloidaldiffusio-phoresis induced by solute gradients [6], or even in the context of biological systems forthe motility of microorganisms [22]. It is therefore necessary to quantify these flows as a functionof the density gradients that generate them for predicting their possible role on other species in thecase of complex fluid mixtures.In the present work, we study in depth the configuration presented in Fig. 1 in a microfluidiccontext, i.e. height H in the (cid:39) µ m range. This experimental configuration has been im-plemented many times in microfluidic devices, either using valves or sliding walls for instance,for various applications such as protein crystallization or biochemical assays [23, 24]. Our maingoal is to quantitatively delineate the range of Rayleigh numbers for which mixing is impacted bybuoyancy in such a configuration, and to predict the laws W vs. T . Configurations similar to thatshown in Fig. 1 have been studied for chemical or civil engineering applications and environmentalissues that involve length scales H ranging typically from 0 . sequestration. Surprisingly, we are not aware of any work that has studied the microfluidic casewhere diffusion cannot be neglected, which also motivated this work. Since our work is related tomicrofluidic applications, we have explored Rayleigh numbers up to Ra = . These high valuesare at the limit of most microfluidic dimensions, but can easily be obtained as soon as H reachesthe millimeter scale, even for dilute solutions. For example, in the numerical application givenpreviously, water and salty water at 1 M, Ra = for H (cid:39) µ m. As we subsequently considerliquid mixtures only, the smallest Schmidt number we explored is Sc = . Such small values can4e observed in the case of the diffusion of small molecules in a low-viscosity solvent, e.g. waterin acetone [31].The present paper is organized as follows. In Sec. II, we present the set of equations modelingthe transport of the solute in the configuration shown in Fig. 1, as well as details about the nu-merical resolutions. In Sec. III, we study the case of a 2D slit for the sake of simplicity, i.e. twoinfinite and parallel plates separated by a thickness H . The numerical data show a rich temporalsuccession of different regimes of solute spreading, that can be captured using analytical asymp-totic solutions, and the analogy with the case of a 2D porous layer [30] is discussed. We finallyaddress in Sec. IV the case of a 3D microfluidic channel for which transverse flows also exist. Wefinally conclude our work in Sec. V and insist on its possible implications. II. MODEL AND NUMERICAL RESOLUTIONA. Model and dimensionless variables
We consider the situation described in Fig. 1(a): a straight microfluidic channel of infinitelength and rectangular cross-section initially filled with a solution at concentration Φ i for Z > Z <
0. For the sake of simplicity, we consider that the kinematic viscosity ν and the interdiffusion coefficient D are constant, and that the density of the solution evolveslinearly with the volume fraction in solute, Eq. (1). Because of this linearity, the problem describedhere can trivially also apply to the interdiffusion between two solutions of different concentrations.Assuming an isothermia of the system and the Boussinesq approximation, the equations gov-erning the solute transport and the velocity field U are: ρ (cid:18) ∂ U ∂ T + U . ∇ U (cid:19) = ρ ν ∆ U − ∇ P + ( ρ ( Φ ) − ρ ) g , (5) ∇ . U = , (6) ∂ Φ ∂ T + U . ∇Φ = D ∆Φ , (7)where P is the pressure deviation from the hydrostatic pressure field for Φ =
0. Boundary con-ditions at the solid walls are the no-slip and the impermeability conditions, U = n . ∇Φ = U ( Z → ± ∞ ) = ∂ Z Φ ( Z → ± ∞ ) = U = Φ = Φ i for Z > Φ = Z < x = X / H , y = Y / H , z = Z / H , t = DT / H , γ = L / H , (8) u = ( H / D ) U , p = H / ( ρ ν D ) P , ϕ = Φ / Φ i . (9)With such definitions, the model given by Eqs. (5-7) reads now:1Sc (cid:18) ∂ u ∂ t + u . ∇ u (cid:19) = ∆ u − ∇ p − Ra ϕ e x , (10) ∇ . u = , (11) ∂ ϕ∂ t + u . ∇ ϕ = ∆ ϕ , (12)where e x is the unit vector along x . With these variables, the initial conditions are: u ( x , y , z , t = ) = and ϕ ( x , y , z , t = ) = H ( z ) , (13)where H ( z ) is the Heaviside function. Boundary conditions are given by: u = n . ∇ ϕ = u ( z → ± ∞ ) = ∂ z ϕ ( z → ± ∞ ) = . (15)To estimate the role of buoyancy on the solute spreading, we first define the cross-sectionaveraged concentration profile by: ϕ ( z , t ) = < ϕ > = γ (cid:90) γ / − γ / (cid:90) ϕ ( x , y , z , t ) d x d y , (16)and the extent of the mixing zone by: w ( t ) = (cid:115) (cid:90) ∞ − ∞ z ∂ ϕ ∂ z d z . (17)In the case of a neutrally-buoyant solute, Ra = u = , (18) ϕ ( x , y , z , t ) = ϕ ( z , t ) = (cid:20) + Erf (cid:18) z √ t (cid:19)(cid:21) . (19)In this case, the width of the interdiffusion zone is given by w = √ t , i.e. W = √ DT with realunits. This is the classical square-root spreading of the solute due to molecular diffusion. ForRa >
0, any deviation from this simple law is a priori a signature of buoyancy-induced dispersion,see Fig. 1(b). 6 . Numerical resolution
Eqs. (10-14) have been solved numerically for two distinct geometries: the 2D case of a mi-crofluidic slit ( γ → ∞ ) and the 3D case of a rectangular micro-channel with a square cross-section( γ = z → ± ∞ have been moved to z = ± λ where λ is a finite distance such that λ (cid:29) w .The 2D numerical simulations have been performed with the commercial software ComsolMultiphysics based on finite elements (Galerkin method). Time discretization is based on implicitBackward Differentiation Formulas, with an adaptive time stepping. Spatial discretization wasachieved by a structured mesh of Lagrangian elements, linear for the pressure and quadratic forthe other variables. The mesh convergence has been thoroughly tested by successive refinements.Computations were made on a workstation with 32 Intel Xeon 2 .
10 GHz processors and 250 GBof RAM.The 3D simulations required the use of in-house made software [33, 34], specifically opti-mized for the simulation of free convection in cavities on parallel architectures and based on amultidomain spectral method. Chebyshev collocation is used for spatial discretization of the threedimensions of space. The pressure-flow coupling is ensured by a projection method that forces thevelocity divergence-free condition. Time integration is performed through a second order temporalscheme combining a Backward Differentiation (BDF2) scheme for the linear terms with an AdamsBashforth extrapolation for the convective terms. Domain decomposition along the z -horizontaldirection is carried out by the Schur complement method for parallelization purposes. Each spa-tial domain is a cube of size 1 in dimensionless units. The spatial resolution in the z -direction hasbeen increased at the first moments of the simulation in order to capture the stiff concentrationgradient and the mesh convergence has been checked by observing the decay of the Chebyshevspectral coefficients. Computations were made in a HPC facility, using from 40 to 360 Intel Xeon2.30 GHz processors.For both 2D and 3D geometries, simulations have been divided into several time intervals inorder to adapt the simulation parameters to the temporal evolution of the flow. For each new timeinterval, the length 2 λ of the spatial domain was extended to take into account the increase of themixing zone. 7 II. THE CASE OF A SLIT
In this section, we study the case of a microfluidic slit, and the different fields in Eqs. (10–12)now depend only on the two variables x and z . We performed numerical simulations for a fixedSchmidt number Sc = , and three different Rayleigh numbers Ra = , 10 , and 10 , over awide range of time scales, from t = − to t (cid:39) . To test the role of the Schmidt number, wealso performed numerical resolutions for the same Rayleigh numbers, Sc = , 10 , 10 and 10 ,and time scales ranging from t = − to t (cid:39) × − as Sc only plays a role at early time scales,see below. A. The case Ra = and Sc = We begin with the case Ra = and Sc = . Figure 2(a) displays several 2D concentrationmaps ϕ ( x , z , t ) at the times shown in Fig. 2(b), see also movie M1 corresponding to these data in theESI. These data evidence a combination of solute spreading by diffusion and buoyancy-inducedadvection.Figure 2(b) displays w ( t ) computed from the numerical resolution for the case Ra = andSc = , along with the diffusion law w = √ t . These data clearly show that the mixing zoneevolves according to the diffusion law expected without buoyancy at small time scales t (cid:28) − ,but also at long time scales t (cid:29) . For intermediate times, the width of the mixing zone issignificantly larger, evidencing the role of the buoyancy-driven advection. The w vs. t behaviorcan be rationalized using different regimes, each with a given power law w ∼ t δ shown in Fig. 2(b).These regimes are presented below in detail, along with self-similar asymptotic solutions for theconcentration profiles ϕ ( z , t ) and the corresponding spreading laws w vs. t . B. Early diffusion regime
We first analyse the transport of solute at short time scales. Figure 3(a) displays the height-averaged concentration profile ϕ ( z , t ) [Eq. (16)] for time scales t ≤ . × − evidencing thespreading of the solute. As shown in Fig. 3(b), all the profiles collapse on a single curve whenplotted against the reduced variable z / √ t . This curve is correctly described by Eq. (19) demon-strating that the transport of the solute is dominated by diffusion for these early time scales, i.e.negligible effect of the flow on the solute transport. As a result, the width of the mixing zone com-8 -6 -4 -2 -2 w t ∼ t / ∼ t ∼ t / ∼ t / ∼ t / z t ∼ / Ra ∼ / Ra ∼ ∼ Ra early diffusionearly advectionlate advection1D dispersionlate diffusion w x (a) (b) w = √ t FIG. 2. For Ra = and Sc = : (a) concentration fields ϕ ( x , z , t ) at the times indicated by the corre-sponding symbols in (b), the z -scale is different in each image. (b) Width of the mixing zone w ( t ) . Themagenta dashed line is given by Eq. (27) and corresponds to early advection with w ∼ t . The yellow dashedline is given by Eq. (29) and corresponds to late advection with w ∼ √ t . The green dashed line is computedfrom the numerical resolution of Eq. (31) and includes both the 1D dispersion regime w ∼ t / and latediffusion w ∼ √ t . The blue dashed line is the diffusion law w = √ t . See also movie M1 corresponding tothese data in the ESI. puted from the 2D data using Eq. (17) is correctly fitted by w = √ t for t ≤ . × − as shown inFig. 2(b).However, a flow exists as u = z -component of the pressure gradient due to buoyancy. Figure 4(a) indeed displaysthe flow driven by this difference of density for t (cid:39) . × − . This gravity current correspondsto a recirculating flow developed within the slit on a length scale ∼
1. Figure 4(b) showing themaximal value of the component u z in the plane z = (cid:39) for time scales t ≥ . × − for Sc = . This plot also shows the same data corresponding to several Schmidt numbers,Sc = , 10 , 10 and 10 . After a transient, the velocity reaches the same plateau value for allSchmidt numbers, except for Sc = for which the plateau is not reached.As the transients depend on Sc and thus on the Grashof number Gr, see Eq. (4), these results9 -3 -5 0 500.20.40.60.81 -5 0 500.20.40.60.81 (b)(a) z ϕ ( z , t ) z/ √ t ϕ ( z , t ) FIG. 3. (a) Height-averaged concentration profiles ϕ ( z , t ) vs. z for time scales ranging from t = − to t ≤ . × − (7 curves, Ra = and Sc = ). (b) Same data plotted against the reduced variable z / √ t ,the dashed line is given by Eq. (19). suggest the existence of an inertial regime corresponding to the development of the recirculatinggravity current through the slit. To better highlight this regime, Fig. 4(c) displays the maximal z -component of the velocity field at z = t for several Rayleigh numbers Ra = ,10 , 10 and several Schmidt numbers Sc = , 10 , 10 and 10 . All the transients collapseon a single curve when times are scaled by 1 / Sc and velocities by Ra. This result is recoveredfrom the Navier-Stokes equation Eq. (10) assuming that the non-linear inertial term u . ∇ u doesnot play any role. The start-up of the flow therefore corresponds simply to the diffusion of themomentum through the slit, expected to take place on a time scaling as ∼ / Sc. By estimatingthe time it takes for the maximal velocity u z to reach 90% of its plateau value, the duration of theinertial regime is t (cid:39) . / Sc, see the vertical dashed line in Fig. 4(c). The smallest Schmidtnumber investigated, Sc = , is an exception. In this case, the steady plateau is not observedbecause inertia is still significant after the end of the diffusion regime, i.e. when advection startsaffecting the solute transport. Turning to real units, the duration of the inertial regime is givenby (cid:39) . H / ν , and lasts only a few tens of milliseconds even for H = µ m and low-viscosity solvents ν = × − m /s. This numerical application shows that such a regime cannotbe observed in most microfluidic experiments, as expected, and that only the Rayleigh number,related to the competition between gravity-induced advection and diffusion, governs the transportof the solute.We now turn to the part of the diffusive regime characterized by a steady gravity current, i.e.10 -6 -4 (b)(a) z t m a x [ u z ( x , z = ) ] -5 -5 -4 -3 -2 (c) m a x [ u z ( x , z = ) ] / R a Sc t x FIG. 4. (a) Velocity vector field superimposed with the concentration field at t = . × − (Ra = ,Sc = ). (b) Maximal value of the component u z in the plane z = t for Ra = and Sc = ( (cid:63) ),10 ( (cid:3) ), 10 ( ◦ ), 10 ( (cid:79) ). The horizontal dashed line is (cid:39) u z ( x , z = ) rescaledby Ra vs. Sc t for Ra = (blue), 10 (magenta), and 10 (cyan) and Sc = ( (cid:63) ), 10 ( (cid:3) ), 10 ( ◦ ), 10 ( (cid:79) ). The horizontal dashed line is estimated using Eq. (22) and given by (cid:39) . t = . after the end of the inertial transient, see the plateau in Fig. 4(b). For the time scales of this regime,neither diffusion nor advection have significantly widened or distorted the concentration field, andit remains close to the initial condition Eq. (13), as evidenced by Fig. 4(a) for the case Ra = .The velocity field after the inertial transient is therefore expected to be the solution of the steadyStokes equation: 0 = ∆ u − ∇ p − Ra H ( z ) e x , (20)where ∇ . u = H ( z ) is the Heaviside function. These equations can be solved analytically,see Appendix A, leading to the two following expressions for the components u x in the plane x = / u z in the plane z = u px ( x = / , z ) = π (cid:90) ∞ sinh (cid:0) k (cid:1) (cid:2) k − (cid:0) k (cid:1)(cid:3) sin ( kz ) k [ k + sinh ( k )] d k , (21) u pz ( x , z = ) = − Ra π (cid:90) ∞ ( x − ) sinh ( kx ) + x sinh ( k − kx ) k [ k + sinh ( k )] d k . (22)These two expressions correctly approximate the velocity profiles in the plateau regime, see theblack lines in Figs. 5(a) and 5(b) for the case Ra = and Sc = . The maximal velocities are11bout max [ u pz ( x , z = )] (cid:39) . x (cid:39) .
81, the plateau value in Fig. 4(c)], and max [ u px ( x = / , z )] (cid:39) . z (cid:39) − . -1 0 1-600-400-2000200400600 0 0.5 1-1000-50005001000 (b)(a) z x u x ( x = / , z ) u z ( x , z = ) FIG. 5. (a) u x ( x = / , z ) vs. z and (b) u z ( x , z = ) vs. x for time scales ranging from t = . × − to 10 − (10 curves, Ra = and Sc = ). The thin black lines are given by Eq. (21) in (a) and by Eq. (22) in (b). C. Early advection
The initial diffusion regime described in Sec. III B ceases when the effect of advection on thethe solute transport becomes non negligible. This transition time corresponds to the departure fromthe diffusion law w = √ t at t (cid:39) × − for the case Ra = and Sc = , see Fig. 2(b). A typ-ical concentration field in this new regime is shown in Fig. 6(a) and has two main characteristics,which can be considered as a definition of the early advection regime: a. The 2D deformation of the concentration field, obviously due to advection, is muchlarger than the diffusive spreading. For this reason, we can neglect diffusion against advection andapproximate the concentration field by an Heaviside function: ϕ ( x , z , t ) = H [ z − z f ( x , t )] , (23)where z f ( x , t ) is the position of a front separating two regions, one where ϕ (cid:39) ϕ (cid:39)
1. In this advection regime, the front motion along z during a small time interval d t reads: d z f ( x , t ) = u z ( x , z f ( t ) , t ) d t . (24)12 -0.1 0 0.1 -0.2 -0.1 0 0.1 0.200.20.40.60.81 -1 -0.5 0 0.5 100.20.40.60.81 x z z (a) (b) ϕ ( z , t ) z/ (Ra t ) (c) z FIG. 6. Concentration fields ϕ at times (a) t (cid:39) . × − and (b) t (cid:39) . × − for Ra = and Sc = .The magenta dashed lines are the deformations estimated by Eq. (25). (c) Average concentration profilesagainst z / ( Ra t ) for time scales ranging from t (cid:39) × − to 1 . × − (9 curves), the dashed line is givenby Eq. (26). Inset: same profiles against z . b. The 2D deformation due to advection is much lower than 1, i.e. than the channel height.The concentration field is thus very close to the initial one, and the velocity field is still givenby Eqs. (21) and (22) corresponding to the plateau observed in Fig. 4(b). We deduce from theseassumptions and Eq. (24) that the front profile can be approximated by: z f ( x , t ) (cid:39) u pz ( x , z = ) t , (25)where u pz ( x , z = ) is the steady velocity field given by Eq. (22). Eq. (25) is consistent with theconcentration field obtained from numerical simulations, as shown in Fig. 6(a).Eqs. (22), (23) and (25) allow the estimation of the height-averaged concentration profile inthis regime: ϕ ( z , t ) = (cid:90) H [ z − u pz ( x , z = ) t ] d x . (26)One can easily show that the above relation is a self-similar function of the variable z / ( Ra t ) .Figure 6(c) compares the theoretical relation Eq. (26) with the numerical simulations, evidencinga reasonably good collapse of the data on the theoretical master curve. Furthermore, the width ofthe mixing zone w ( t ) defined by Eq. (17) can be estimated using Eq. (26), leading to: w (cid:39) . t . (27)This behavior is plotted in Fig. 2(b) and correctly fits the data obtained from the numerical resolu-tion of the 2D model from t = . × − to 2 × − .13t later time, significant discrepancies are observed between the theoretical formula Eq. (25)and the numerical simulation, see for instance Fig. 6(b) for a comparison at t = . × − . Indeed,the deformation of the concentration field is of the order of 1 for such time scales and Eq. (22) canno longer be used for the estimation of the velocity field. It marks the end of the early advectionregime. D. Late advection
At later time scales, solute spreading by diffusion still remains negligible as compared to soluteadvection, but the 2D deformation of the concentration fields is now much larger than 1, see forinstance Fig. 7(a) showing a snapshot at t = . × − . These data along with the flow field shown -0.3 -0.2 -0.1 0 0.1 0.200.20.40.60.81 -3 0 3 -1000100 -6 -4 -2 0 2 4 600.51 -50050 x zxx (a) (d) (b)(c) ϕ ( z , t ) ξ = z/ (cid:31) ˜ Dt z ϕu z u x FIG. 7. (a) Concentration field ϕ , and components (b) u z , (c) u x at t = . × − for Ra = and Sc = .The magenta dashed line is given by the solution of Eq. (28) in the case of a slit. (c) Average concentrationprofiles ϕ ( z , t ) vs. ξ = z / √ ˜ Dt for time scales ranging from t (cid:39) .
003 to 0 .
02 (9 curves), the dashed line is ψ ( ξ ) solution of Eq. (28). Inset: same profiles against z . in Figs. 7(b) and 7(c), evidence the reciprocal exchange of the solution and the solvent separated bya diffuse pseudo-interface. For fully negligible diffusion, this regime commonly referred to as theviscous lock-exchange problem has been widely studied in the literature. In such a configuration,many groups predicted that the extent of the spreading of the two fluids scales as w ∼ √ ˜ Dt where ˜ D is an effective diffusion coefficient. The square-root behavior arises from the competition betweenbuoyant forces ( ∼ ρ β Φ i gH / W ) and viscous forces ( ∼ ρ ν ˙ W / H ), leading to the dimensionless14caling law w ∼ Ra t [28]. The effective diffusion coefficient ˜ D therefore only depends on Raand the geometry, and it has been calculated in various cases: porous medium, circular tube, butalso rectangular channel and slit [27–30]. The main idea of these works is to compute the shape ψ ( z , t ) of the pseudo-interface separating the two fluids, assuming large deformations and thus aquasi-parallel flow along the z -axis (lubrication approximation). With such approximations, onecan show that ψ ( z , t ) admits a self-similar shape ψ ( ξ ) with ξ = z / √ ˜ Dt , solution of: − ξ d ψ d ξ = ξ (cid:18) f ( ψ ) d ψ d ξ (cid:19) , (28)with f ( ψ ) = ψ ( − ψ ) and ˜ D = Ra / et al. in Ref. [27]. Figures 7(a–c) show thissolution for the corresponding time t = . × − superimposed with both the concentration fieldand the velocity field. These data show a reasonable agreement, confirming the negligible impactof diffusion on the solute transport.To better describe this regime, Fig. 7(d) displays the average profiles ϕ ( z , t ) obtained from the2D model at Ra = and Sc = for time scales ranging from t (cid:39) .
003 to 0 .
02. This plotshows that all the data almost collapse on a single curve when plotted against ξ = z / √ ˜ Dt , whichis correctly described by ψ ( ξ ) solution of Eq. (28). In this regime, one can again compute thewidth of the mixing zone defined by Eq. (17) using ψ ( ξ ) , leading to: w (cid:39) . √ Ra t . (29)This square-root spreading is plotted in Fig. 2(b) and accounts well for the numerical data w ( t ) obtained for Ra = and time scales ranging from t (cid:39) − to (cid:39) − . E. 1D dispersion and late diffusion
For time scales t ≥ O ( ) , diffusion almost homogenizes the solute over the height of the slit,and the transport of the solute cannot be described by only advection, as revealed by Fig. 8(a)showing the concentration field at t = .
5. In this regime, the transport is fully controlled bythe coupling between solute diffusion along the channel height and buoyancy-driven advectionalong the channel main axis. This regime has already been described in the literature since thepioneering work of Chatwin and Erdogan, who studied the Taylor-Aris dispersion of a buoyantsolute in a pressure-driven flow [35], see also [36–40] and the review of Young and Jones on shear15 -100 0 100 -5 0 500.51 -1500 0 1500 -20020 -30 -20 -10 0 10 20 3000.51 -101 x z xx (a) (d) (b)(c) ϕ ( z , t ) ϕ ( z , t ) (e) ( α / / Ra / ) z/t / z/ √ t zz FIG. 8. (a) Concentration field ϕ , (b) component u z and (c) component u x at t = . = and Sc = ).The black lines in (a) are isoconcentration lines. (d) Average concentration profiles ϕ ( z , t ) in the dispersionregime plotted against η given by Eq. (34), t ranges from t = . × (10 curves). The dashed lineis Eq. (33). Inset: same profiles against z . (e) Rescaled concentration profiles ϕ ( z , t ) in the late diffusionregime, t = up to 1 . × (10 curves). The dashed line is Eq. (19). Inset: same profiles against z . dispersion [41]. In this regime, the extent w of the concentration gradient along z is large ( w (cid:29) u x (cid:28) u z , see Figs. 8(b) and 8(c)], and the variations ofthe concentration along x are small [see the isoconcentration lines in Fig. 8(a)]. One can thereforeuse the lubrication approximation to show that the density gradient along z drives a flow following: u z ( x , z , t ) = − Ra12 ∂ ϕ ∂ z x ( x − )( x − ) , (30)see Appendix B. This flow adds a contribution to the dispersion which scales as ∼ u z , as for theclassical Taylor-Aris dispersion in a Poiseuille flow. More rigorously, one can demonstrate thatthe average concentration ϕ ( z , t ) obeys the 1D dispersion equation: ∂ ϕ ∂ t = ∂∂ z (cid:18) D eff ∂ ϕ ∂ z (cid:19) , (31)with: D eff = + α (cid:18) Ra ∂ ϕ ∂ z (cid:19) , (32)and α = t (cid:63) = ( α / Ra ) t and z (cid:63) =( √ α / Ra ) z . It is then solved numerically with the initial condition ϕ ( z , t = ) = H ( z ) to computethe solution for any Rayleigh number and α value. Figure 2(b) shows that the width of the mixingzone w ( t ) defined by Eq. (17) and computed from the numerical resolution of Eq. (31), perfectlymatches the data obtained from the full 2D model for time scales t ≥ .
1. The component u z at z = ϕ ( z , t = ) = H ( z ) has been studied by Macleanand Alboussi`ere [39] who provided asymptotic approximations of the solution. When buoyancydominates the transport of the solute, i.e. D eff (cid:29)
1, Eq. (31) admits the self-similar solution: ϕ = + η (cid:115) √ π − η + π arcsin (cid:18) √ πη × / (cid:19) , (33)with η given by: η = α / √ Ra zt / . (34)Eq. (33) is valid for η ≤ / ( π √ ) [39]. Figure 8(d) displays this asymptotic self-similar solutionfor Ra = along with the data computed from the 2D model, evidencing a very good agreement.In this regime, one can compute the width of the mixing zone leading to: w (cid:39) (cid:114) Ra2 π (cid:18) t α (cid:19) / , (35)thus following w ∼ t / as shown in Fig. 2(b).At later times, the density gradient continuously decreases as the solute is continuously dis-persed along the channel, and diffusion dominates again the transport of the solute, i.e. D eff (cid:39) w = √ t as shown in Fig. 2(b), despite the buoyancy-driven flow along the slitstill given by Eq. (30). As explained in Introduction, even if this flow has no effect on the solutegradient that generates it, it still exists and may have an effect on less mobile species in the caseof complex liquid mixtures. 17 - early diffusion w (cid:39) √ t ¯ u z ( z = , t ) (cid:39) . t I → II (cid:39) / Ra II- early advection w (cid:39) . t ¯ u z ( z = , t ) (cid:39) . t II → III (cid:39) / RaIII- late advection w (cid:39) . √ Ra t ¯ u z ( z = , t ) (cid:39) . (cid:113) Ra t t III → IV (cid:39) . w (cid:39) (cid:113) Ra2 π (cid:0) t α (cid:1) / ¯ u z ( z = , t ) (cid:39) √ Ra192 √ π (cid:0) α t (cid:1) / t IV → V (cid:39) ( π ) α Ra V- late diffusion w (cid:39) √ t ¯ u z ( z = , t ) (cid:39) √ π Ra √ t TABLE I. Transport regimes and transition times for a slit, α = F. Dispersion diagram
In the previous paragraphs, we identified a sequence of regimes of transport of the solute, whichare summarized in Table I. For each regime, one can compute the typical longitudinal velocity at z = u z ( z = , t ) = − (cid:90) / u z ( x , z = , t ) d x , (36)using in particular Eq. (22) for Regimes I and II, Ref. [27] for Regime III, and Eq. (30) alongwith Eq. (33) [resp. Eq. (19)] for Regime IV (resp. Regime V). Results are displayed in the thirdcolumn of Table I.The transition times between these regimes are estimated by matching the different spreadinglaws w vs. t leading to the values provided in Table I. These definitions lead to different scalinglaws with Ra which are also reported in Fig. 2. These transition times allow us to construct thediagram presented in Fig. 9(a) in the plane Ra vs. t . The vertical dashed-dotted lines in Fig. 9(a)corresponding to the extent of the inertial regime show that the early advection regime II is alsocoupled for large Ra and small Sc to the momentum diffusion across the slit, see Sec. III B and inparticular Fig. 4(b).This diagram also shows that the effect of buoyancy vanishes for Ra (cid:46) . To illustrate thispoint, Fig. 9(b) reports the spreading laws w vs. t computed from the numerical resolution of thefull 2D model for several Rayleigh numbers Ra = , 10 , 10 and Sc = . These data clearlyreveal that buoyancy has little effect on the transport of the solute at all time scales for Ra = .More quantitatively, the ratio w / √ t reaches a maximum of only (cid:39) .
24 at t (cid:39) . = .This result answers the question initially asked in Introduction as it allows to assess a numericalvalue to Ra corresponding to negligible buoyancy in a microfluidic slit at all time scales, at leastregarding the active solute that generates the gravity current.18 -6 -4 -2 diffusion D d i s p e r s i o n e a r l y a d v e c t i o n l a t e a d v e c t i o n t R a S c = S c = e a r l y d i ff u s i o n l a t e d i ff u s i o n w = √ t -6 -4 -2 -2 w t (a) (b) FIG. 9. (a) Diagram of the different regimes in the plane Ra vs. t . The transition times t vs. Ra are given inTable I. The vertical dotted lines correspond to the extent of the inertial regime t (cid:39) . / Sc for Sc = and Sc = , see Sec. III B. (b) Width of the mixing zone w ( t ) for Sc = and Ra = (red), 10 (blue), 10 (magenta). The black dashed line is the purely diffusive spreading w = √ t . The symbols are thetransition times given in Table I for Ra = and Ra = . The regimes of early diffusion and early advection are only visible for high Rayleigh numbers,Ra ≥ , and small time scales, t < − . For most microfluidic configurations investigatingmolecular solutes, these regimes might be difficult to observe even in a thick slit, as t = − in theabove diagram does not exceed a few seconds for H = µ m and D ≥ − m /s. On the otherhand, both regimes of late advection and 1D dispersion should be easily observed as the transitiontime between these two regimes t III → IV is a few minutes for D (cid:39) − m /s and H = µ m.In the case of colloidal dispersions, the regimes of early diffusion and early advection might beobservable as much lower D values lead to much longer time scales. For instance, t = − is nowa few minutes for H = µ m and D = − m /s corresponding to colloids of radius 100 nmdispersed in water. G. Analogy with the case of a 2D porous layer
Szulczewski and Juanes [30] studied in a different context (geological sequestration of CO inan aquifer), a problem similar to the one described in Fig. 1, but considering a vertically confinedporous layer of thickness H . Their theoretical model is also based on Eq. (7) to describe the19olute transport and Eq. (6) for the overall mass conservation, but the pore velocity field U followsDarcy’s law given by: U = − κρ ν ε [ ∇ P − ( ρ ( Φ ) − ρ ) g ] , (37)where κ is the permeability of the permeable rock and ε its porosity (see Sec. II A for the othernotations). Notice that Navier-Stokes equations Eq. (5) turns to Stokes equations when inertiais neglected, that does not reduce to Darcy’s law Eq. (37). Indeed, Stokes equations include thediffusion of the momentum over the scale H of the microfluidic channel, that brings a fundamentaldifference with Darcy’s law.The characteritic velocity resulting from the Darcy’s law Eq. (37) reads U D = β Φ i g κν ε , (38)and the corresponding Rayleigh number comparing diffusion and advection by the gravity currentis: (cid:102) Ra = β Φ i g κ H ν ε D ∼ U D HD , (39)thus highlighting a scaling law with the thickness different from the microfluidic case, (cid:102) Ra ∝ H vs. Ra ∝ H . Interestingly, the range of Rayleigh numbers (cid:102) Ra involved in the context of CO sequestration [30] still corresponds to the range of Ra studied in the present work.Despite the differences between these two models, Szulczewski and Juanes also reported fivedistinct transport regimes in the porous rock that have strong similarities to those reported inFig. 2(a) for a microfluidic slit, see their names in Table II. In particular, the concentration profilesdepend on the same self-similar variables in each regime, and the different transition times (dis-played in the last column of Table I for our model) obey the same scaling laws with the Rayleighnumber. We believe that these similarities are related to the linearity of the Darcy and Stokes equa-tions in both problems. Nevertheless, since the velocity fields are different in both configurations,we expect possibly different prefactors for the scaling laws. To confirm this point, we estimate thesolute flux across the interface z = f ( t ) = − (cid:90) (cid:20) ϕ u z − (cid:18) ∂ ϕ∂ z (cid:19)(cid:21) x , z = d x , (40)as Szulczewski and Juanes also computed this quantity in each regime [30]. f is estimated usingEqs. (19) in Regimes I and V, Eq. (22) in Regime II, following Ref. [27] in Regime III, and using20 icrofluidic slit, this work 2D porous layer [30]I- early diffusion f (cid:39) ( π t ) / I- early diffusion ˜ f (cid:39) ( π t ) / II- early advection f (cid:39) . f (cid:39) . (cid:102) RaIII- late advection f (cid:39) . Ra / t / III- straight-line slumping ˜ f (cid:39) . (cid:102) Ra / t / IV- 1D dispersion f (cid:39) . Ra / t / IV- Taylor slumping ˜ f (cid:39) . (cid:102) Ra / t / V- late diffusion f (cid:39) ( π t ) / V- late diffusion ˜ f (cid:39) ( π t ) / TABLE II. Solute flux defined by Eq. (40) for a microfluidic slit and in a 2D porous layer [30].
Eq. (31) in Regime IV. As shown in Table II, f for the microfluidic slit and ˜ f for the porouslayer show the same scaling laws with t and Ra (resp. (cid:102) Ra). With the exception of the diffusionRegimes I and V, the large differences between the numerical prefactors (up to two orders ofmagnitude) confirm that both problems are fundamentally different. Relevant prefactors must beconsidered in potential comparisons with experimental works.
IV. THE CASE OF A 3D MICROFLUIDIC CHANNEL
We now consider the case of a microfluidic channel with a square cross-section, i.e. γ = L / H = = , 10 , and10 with a fixed Schmidt number Sc = . A. Transport regimes in the 3D case
Movie M2 supplied in the ESI shows the concentration and velocity fields in the planes y = z = = and Sc = .This movie helps to identify the same succession of regimes as for the slit case. From the fullnumerical data, we computed again the width of the mixing zone w ( t ) using Eqs. (16) and (17) forthree different Rayleigh numbers Ra = , 10 and 10 and the same Schmidt number Sc = ,see Fig. 10. The dispersion curves obtained from the 3D model follow the same trends as the 2Dones shown in Fig. 2(b). These data lead to the same conclusion as for the slit case: buoyancyhardly affects the solute transport at all time scales for Ra ≤ (the maximal value of the ratio w / √ t is only (cid:39) .
14 at t (cid:39) . = ).These observations evidence that the description of the regimes presented above applies again21 -4 -2 -2 w t ∼ t / ∼ t / ∼ t / w = √ t FIG. 10. Width of the mixing zone w ( t ) obtained in the case of a channel with a square cross-sectionfor Sc = and Ra = (red), 10 (blue), 10 (magenta). Asymptotic models for Ra = : the yellowdashed line given by Eq. (41) corresponds to the regime of late advection; the green dashed line is computedfrom the numerical resolution of Eqs. (31) and (32) with α (cid:39) w = √ t . for the 3D case. Thereafter, we will not re-describe the regimes of early diffusion and earlyadvection, in particular because they are short or even hardly observable in most microfluidicexperimental configurations. With respect to the late advection regime, the analysis reported inSec. III D can easily be adapted to the 3D case. The shape of the pseudo-interface separating thesolution and the solvent is computed using Eq. (28) with f ( ψ ) and ˜ D corresponding to a channelwith a square cross-section, see Ref. [27] for details. We then calculated the width w ( t ) from thetheoretical profiles, leading to: w (cid:39) . √ Ra t . (41)This prediction fits well the data obtained from the 3D numerical simulation at Ra = for timescales ranging from t (cid:39) × − to (cid:39) × − , see Fig. 10. The numerical prefactor in Eq. (41)(0 . . t ≥ O ( ) , Movie M2 evidences that diffusion almost homogenizes theconcentration over the cross-section of the channel. This is illustrated by Figs. 11(a) and 11(d)22howing the concentration field in the planes y = z = t (cid:39) .
7. Figure 11, displayingalso the components of the velocity field in the same planes and at the same time, evidences aquasi-parallel flow along z (due to the large extent w of the longitudinal density gradient), and asecondary transverse flow as revealed by the components u x and u y in the z = -10010 -50 0 5000.51 -0.500.5 -0.5 0 0.500.51 -0.5 0 0.500.51 -0.500.5 -0.5 0 0.5 -0.500.5 (a)(b) x z xx xx y y (d)(f) (e)(g) (c) FIG. 11. 3D case, Ra = , Sc = , and t (cid:39) .
7. (a) Concentration field ϕ ( x , y = , z , t ) . (b) Velocityfield u z ( x , y = , z , t ) and (c) u x ( x , y = , z , t ) . (d) ϕ ( x , y , z = , t ) , (e) u z ( x , y , z = , t ) , (f) u x ( x , y , z = , t ) , and(g) u y ( x , y , z = , t ) . The black lines in (a) and (d) are isoconcentration lines. All these results suggest, as for the 2D case of the slit, the existence of a dispersion regimedescribed by a 1D model, see Eq. (31). However, the 3D case deserves particular attention because(i) to our knowledge, the expression Eq. (32) of the dispersion coefficient cannot be found in theliterature for a rectangular channel. It has been already calculated for a circular tube [35] and aslit [41], and we derive it for a rectangular cross-section below in Sec. IV B; (ii) As mentioned byChatwin and Ergogan [35] and evidenced in Fig. 11, a transverse flow exists in a 3D geometry,which might hinder the validity of the 1D dispersion model. This issue is addressed below inSec. IV C.
B. 1D dispersion regime in the 3D case
The fact that w (cid:29) ϕ ( x , y , z , t ) = ϕ ( z , t ) + ϕ ( x , y , z , t ) , (42)with ϕ given by Eq. (16) and ϕ (cid:28) ϕ . Averaging the transport equation Eq. (12) over the cross-section of the channel leads to: ∂ ϕ ∂ t + ∂ < u z ϕ > ∂ z = ∂ ϕ ∂ z , (43)similarly to Eq. (B2) for the case of a slit. Subtracting this relation to Eq. (12) gives: ∂ ϕ ∂ t + u z ∂ ϕ ∂ z + u . ∇ ϕ − ∂ < u z ϕ > ∂ z = ∂ ϕ ∂ x + ∂ ϕ ∂ y + ∂ ϕ ∂ z , (44)as Eq. (B3). For t (cid:29) w (cid:29)
1, Eq. (44) yields at leading order: u z ∂ ϕ ∂ z + u y ∂ ϕ ∂ y + u x ∂ ϕ ∂ x (cid:39) ∂ ϕ ∂ x + ∂ ϕ ∂ y . (45)In the 2D case of a slit, u y = u x ∼ u z / w (cid:28) u z . Theterm u x ∂ x ϕ in the above relation is thus negligible, and the derivation of the dispersion equationEq. (31) is straightforward, see Appendix B. In the 3D case nevertheless, the continuity equationdoes not make it possible to relate the scales of u y and u x to u z , as there can exist (possibly large)secondary transverse flows verifying ∂ x u x + ∂ y u y =
0. It is therefore a priori not possible to derivethe dispersion equation Eq. (31) unless one neglects the terms u x ∂ x ϕ and u y ∂ y ϕ in Eq. (44). Thispoint had been mentioned by Chatwin and Erdogan who studied the classical Taylor-Aris problemof the dispersion of a buoyant solute flowing in a circular tube [35]. They even showed that thelateral mixing of the solute due to these secondary transverse flows could lead to a decrease of theoverall solute dispersion in a pressure-driven flow for some range of the Rayleigh number, see alsoRefs. [36, 37].The data shown in Fig. 11 along with Movie M2 in the ESI evidence that the maximal mag-nitude of these secondary flows is about u x (cid:39) u y (cid:39) = , and we will assume as a firststep that they do not significantly affect the transport of the solute. This assumption allows us toneglect the terms u y ∂ y ϕ and u x ∂ x ϕ in Eq. (44) as compared to the diffusive terms ∂ x ϕ and ∂ y ϕ ,leading to: u z ∂ ϕ ∂ z (cid:39) ∂ ϕ ∂ x + ∂ ϕ ∂ y , (46)similarly to Eq. (B4) for the slit. 24o compute the longitudinal flow u z , we first assume that it is described by the lubricationapproximation because the extent of the mixing zone is large ( w (cid:29) ∂ u z ∂ x + ∂ u z ∂ y = ∂ p (cid:96) ∂ z , (47)0 = ∂ p (cid:96) ∂ y , (48)0 = ∂ p (cid:96) ∂ x + Ra ϕ , (49)where p (cid:96) is the pressure field associated to this longitudinal flow. After integrating Eqs. (48)and (49) and inserting the resulting pressure field in Eq. (47), one finds: ∂ u z ∂ x + ∂ u z ∂ y = Ra ∂ ϕ ∂ z (cid:18) − x (cid:19) , (50)the term 1 / < u z > = u z and ϕ , and finally to com-pute the dispersion term < u z ϕ > in Eq. (43) leading ultimately to the same dispersion equationEq. (31) but with a different α value for D eff in Eq. (32). Appendix C reports solutions of Eqs. (46)and (50) using Fourier series for a rectangular channel of arbitrary aspect ratio γ , see Eqs. (C2)and (C3). These calculations allow us to compute the numerical values of α as a function of γ ,see Fig. 15 in Appendix C. For a square cross-section, α (cid:39) α tends towards the valuederived for the slit α = u z and ϕ in the plane z = t (cid:39) .
7. These data wellmatch the numerical data reported in Fig. 11(d) and (e) at the same time, see also the compar-isons along the line y = z = ϕ ( z , t ) are correctly predicted by the 1D model (data not shown). Theseprofiles make it possible to compute the width of the mixing zone w ( t ) [Eq. (17)] plotted inFig. 10 along with the data obtained from the 3D numerical simulation for Ra = . The1D dispersion approach accounts well for the spreading of the solute for t ≥ .
1. Moreover,Eq. (35) giving the w ∼ t / behavior, allows us to compare the 2D and 3D configurations in thisregime. Eq. (35) indeed indicates that the ratio w ( t ) / w ( t ) between the 2D and the 3D cases is ( α / α ) / (cid:39) ( / ) / (cid:39) . ≤
1, again due to the side walls increasing viscousforces in the 3D case. 25 -0.03-0.02-0.0100.010.020.03 -0.5 0 0.500.20.40.60.81 -15-10-5051015 (a) (b) (c) (d) x y y x ϕ ( x , , ) u z ( x , , ) FIG. 12. (a) u z ( x , y , z = ) and (b) ϕ ( x , y , z = ) calculated using Eqs. (C2) and (C3) at t (cid:39) .
7, see thecorresponding numerical data shown in Fig. 11(d) and (e). The longitudinal density gradient is calculatedby the numerical resolution of the dispersion equation Eq. (31). (c) and (d) show the comparisons alongthe line y = z = C. Transverse flows and validity range of the 1D dispersion equation
The comparisons shown in Sec. IV B evidence the validity of the 1D dispersion equation todescribe the overall solute transport at least for Ra (cid:46) , although the transverse flow is notconsidered in the 1D model. In the following, we derive the expression of the transverse flow withthe assumption that it does not significantly affect the concentration field. Then, we determine thecritical Rayleigh number for which this hypothesis no longer holds.Our calculations make it possible to compute the transverse variations in concentration ϕ , seeFig. 12(b). These variations evidence density gradients along y due to the presence of the lateralwalls which impact the longitudinal flow, and thus the solute distribution. These density gradientsare responsible for a secondary flow because they cause a pressure gradient along y . We computethis flow assuming that it is locally invariant along z because of the large extent of the mixing zone( w (cid:29)
1) and solenoidal, i.e. ∂ y u y + ∂ x u x =
0, to ensure the global mass conservation. This flow isthus solution of: 0 = ∂ p t ∂ z , (51) ∂ u y ∂ x + ∂ u y ∂ y = ∂ p t ∂ y , (52) ∂ u x ∂ x + ∂ u x ∂ y = ∂ p t ∂ x + Ra ϕ , (53)26here p t is the pressure field associated to this 2D solenoidal flow. As in the work of Chatwin andErdogan [35], the global flow is therefore assumed to be the superposition of the longitudinal flowdescribed by Eqs. (47-49) with the secondary transverse flow determined by Eqs. (51-53). -0.5 0 0.500.20.40.60.81 -0.5 0 0.500.20.40.60.81 -10 0 10 -0.5 0 0.500.20.40.60.81 -10 0 10 x (a) (b) (c) x y y y FIG. 13. Components u x (a) and u y (b) scaled by 10 − (cid:16) Ra ∂ϕ ∂ z (cid:17) computed from the resolution of thebiharmonic equation Eq (D2). (c) Corresponding transverse flow field in the plane x - y . This solenoidal flowfield results from the density gradients given by ϕ , see Eq. (D2). An analytical approximation of the solution of Eqs. (51-53) is computed in Appendix D for arectangular cross-section of arbitrary aspect ratio γ using a method described by Shankar et al. [44, 45]. Figure 13 shows the theoretical prediction of this secondary transverse flow for γ = u x and u y are:max ( u x ) (cid:39) max ( u y ) (cid:39) × − (cid:18) Ra ∂ ϕ ∂ z (cid:19) . (54)To go a step further into the comparison, Fig. 14 displays the maximal values of the compo-nents u x , u y and u z in the transverse plane z = t obtained from the 3D model for threeRayleigh numbers, Ra = , 10 , and 10 . The theoretical prediction given by Eq. (C2) for thelongitudinal flow u z correctly fits the numerical data even at Ra = for t ≥
1. The theoreticaltransverse components given by Eq. (54) accounts well for the data at Ra = and Ra = [with the gradient of ϕ at z = = . For this Rayleigh number, u x (cid:39) u y (cid:39) -4 -2 -4 -2 -4 -2 -2 -4 -2 (a) (c) (b) t t t m a x [ u x , y , z ( x , y , z = ) ] FIG. 14. Maximal value of the components u z (red), u x (blue) and u y (magenta) in the plane z = = , (b) 10 , and (c) 10 (Sc = for these three cases). Theblack lines are the prediction given by Eq. (C2), and the dotted line the prediction given by Eq. (54). Inthese theoretical predictions, the longitudinal density gradient ∂ z ϕ at z =
1D dispersion regime ( t ≥ z = ( u x ) (cid:39) max ( u y ) (cid:39) . × − Ra √ t , (55)for a square cross-section ( α (cid:39) ( u x ) (cid:39) max ( u y ) ≤ t ≥ . × − Ra , (56)about t ≥ . = , in a remarkable agreement with the data of Fig. 14(c) showing adiscrepancy between the numerical solutions ( u x , u y ) and the predictions given by Eq. (54) for t ≤
10. Nevertheless, the transverse flow does not yet significantly change the longitudinal dispersioneven for Ra = at t > mix the solute laterally, as also shown by Chatwin and Erdogan in adifferent context, the Taylor-Aris dispersion of a buoyant solute in a pressure-driven flow [35–37].In most microfluidic experimental configurations, these effects can a priori be neglected becauseRa ≤ . V. CONCLUSIONS
In the present work, we have studied in detail the impact of buoyancy on solute spreading intwo distinct microfluidic geometries: a 2D slit and a microchannel with a square cross-section,in particular through analytical predictions fully validated by precise numerical resolutions of thetransport equations. One of the main results of our study is to show that for Ra ≤ , solutal freeconvection does not impact solute diffusion at all time scales. Beyond this result, our theoreticalpredictions give also for larger Rayleigh numbers, the time scales (or density gradients) for whichbuoyancy no longer impacts molecular diffusion, see the diagram in Fig. 9. Moreover, these sametheoretical predictions allow to estimate analytically the gravity currents, whatever their role onsolute transport. It is worth remembering that these flows may impact the transport of other speciesdispersed in the flow, even though they do not affect the gradients of concentration of the active species that generate them. As an example, for the experimental case mentioned in Introduction,interdiffusion between water and a NaCl aqueous solution at 1 M in a microfluidic slit of height H = µ m, free convection is not expected to impact the diffusive mixing because Ra (cid:39) u z (cid:39) µ m/s in theearly regimes of diffusion and advection, and still ¯ u z (cid:39) µ m/s for the time scale T = H / D (seeTable I), and could significantly advect less mobile species dispersed in the solutions. The precisecontrol of transport conditions in microfluidic geometries thus possibly opens the way to flowcontrol induced by solute gradients.For the 3D case of a rectangular cross-section channel, our work brings for the first time (toour knowledge) estimates of the 1D dispersion coefficient describing the transport of the soluteat long time scales for any aspect ratio γ = L / H . Our work also highlights a subtle point relatedto 3D geometries: the order of magnitude of the transverse flows cannot be determined fromthe lubrication approximation alone. In the case studied here, these flows, induced by transverse29ensity gradients, remain moderate up to Ra = , and the overall solute spreading is correctlydescribed by a 1D dispersion equation. However, our work predicts that these flows could playa role at higher Rayleigh numbers for experimental situations outside the field of application ofmicrofluidics.It could also be relevant to study more in details the case of shallow channels commonly en-countered in microfluidic applications, and in particular to study more finely the transition betweena channel with a large aspect ratio, γ (cid:29)
1, and the slit. Indeed, the case γ (cid:29) ∼ γ ,see Ref. [43] investigating this issue for the case of the Taylor-Aris dispersion.Finally, we considered in our work the case of an ideal binary solution in the framework of theBoussinesq approximation, see Eqs. (5–7). Microfluidic technologies allow a very fine control ofthe transport conditions (especially mass and momentum), and thus to study interdiffusion in morecomplex mixtures. It would then be useful to go beyond the model described by Eqs. (5–7) to in-clude this complexity: change in viscosity and diffusion coefficient as a function of concentration,role of the reference frame (volume velocity / mass velocity) [46, 47], other transport mecha-nisms (e.g. diffusio-osmosis), etc. In this context, we hope that our work will make it possible todisentangle the role played by solutal free convection from other transport phenomena. ACKNOWLEDGEMENTS
This work was performed using HPC resources from the ”M´esocentre” computing center ofCentraleSup´elec and ´Ecole Normale Sup´erieure Paris-Saclay supported by CNRS and R´egionˆIle-de-France (http://mesocentre.centralesupelec.fr/). JBS also thanks Y. Hallez for discussionsconcerning the lock-exchange problem and ANR OSMOCHIP (ANR-18-CE06-0021) as well asSolvay and CNRS for funding.
Appendix A: Buoyancy-driven flow at early stage for the case of a slit
After a transient corresponding to the diffusion of the momentum across the slit, the velocityfield is solution of the Stokes equation Eq. (20). Introducing the stream function ψ ( x , z ) definedby: u x = ∂ ψ∂ z and u z = − ∂ ψ∂ x , (A1)30q. (20) is equivalent to: ∆ ψ = Ra δ ( z ) , (A2)with δ ( z ) the Dirac function and ∆ the biharmonic operator in 2D. The no-slip boundary condi-tions impose: ∂ ψ∂ z ( x = , z ) = ∂ ψ∂ x ( x = , z ) = . (A3)The flow is expected to vanish far from z = ψ ( x , z → ± ∞ ) =
0. Because d ψ = − u z d x + u x d z , integrating d ψ over a contour along x = x = z → ± ∞ imposes: ψ ( x = , z ) = . (A4)We define the Fourier transform of ψ ( x , z ) by:˜ ψ ( x , k ) = √ π (cid:90) ∞ − ∞ d z ψ ( x , z ) e ikz . (A5)Eq. (A2) turns to: ∂ ˜ ψ∂ x − k ∂ ˜ ψ∂ x + k ˜ ψ = Ra √ π . (A6)The solution of the ordinary differential equation is:˜ ψ ( x , k ) = Ra k + k ( x − ) cosh ( kx ) − kx cosh ( k − kx ) + sinh ( k ) − sinh ( kx ) − sinh ( k − kx ) √ π k ( k + sinh ( k )) , (A7)for the above boundary conditions. The stream function is then computed from the inverse Fouriertransform of ˜ ψ ( x , k ) and the velocity field is found using Eqs. (A1). These calculations lead inparticular to Eqs. (21) and (22) given in Sec. III C. Due to the symmetry of the equations along the x = / u x occurs at x = / z values. Appendix B: Derivation of the advection-dispersion equation for the case of a slit
We assume: ϕ ( x , z , t ) = ϕ ( z , t ) + ϕ ( x , z , t ) , (B1)with ϕ given by Eq. (16) and ϕ (cid:28) ϕ . Averaging the transport equation Eq. (12) over the heightof the slit with the help of the continuity relation Eq. (11), leads to: ∂ ϕ ∂ t + ∂ < u z ϕ > ∂ z = ∂ ϕ ∂ z . (B2)31ubtracting this relation to Eq. (12) yields: ∂ ϕ ∂ t + u z ∂ ϕ ∂ z + u . ∇ ϕ − ∂ < u z ϕ > ∂ z = ∂ ϕ ∂ x + ∂ ϕ ∂ z . (B3)If we now assume that the scale w of the gradient along z verifies w (cid:29)
1, the continuity equationimposes u x ∼ u z / w , and Eq. (B3) reduces to: u z ∂ ϕ ∂ z (cid:39) ∂ ϕ ∂ x , (B4)provided that t (cid:29) ϕ (cid:28) ϕ . Similarly, by neglecting the inertial term in Eq. (10) andassuming again w (cid:29)
1, one find Eq. (30) for the component u z , and u x given by the continuityequation Eq. (11). Integration of Eq. (B4) along with the impermeability boundary condition and < ϕ > = ϕ = − Ra1440 (cid:18) ∂ ϕ ∂ z (cid:19) (cid:16) x − x + x − (cid:17) . (B5)Inserting Eq. (B5) into Eq. (B2) leads to the dispersion equation Eq. (31) with the dispersioncoefficient given by Eq. (32). Comparison of the solution of Eq. (31) with the data obtained fromthe full 2D model shows that the 1D dispersion model is actually valid as soon as t ≥ Appendix C: Derivation of the advection-dispersion equation for the case of a microfluidic channelwith a rectangular cross-section
The solution of Eq. (50) is found by a Fourier sum noticing first that:12 − x = ∞ ∑ n = , even n π sin ( n π x ) for x ∈ ] − [ . (C1)The Fourier series representing the solution of Eq. (50) along with the no-slip boundary conditionsat the solid walls is then: u z = Ra ∂ ϕ ∂ z ∞ ∑ n = , even ( n π ) (cid:18) cosh ( π ny ) cosh ( π n γ / ) − (cid:19) sin ( n π x ) , (C2)with γ = L / H , see Ref. [48] for a similar calculation of the pressure-driven velocity profile in arectangular channel. The Fourier series representing ϕ can now be found using Eq. (46) alongwith the impermeability boundary conditions at the solid walls and the constraint < ϕ > = ϕ = π (cid:18) ∂ ϕ ∂ z (cid:19) ∞ ∑ n = , odd ∞ ∑ p = , even cos ( n π x ) p ( p − n ) (cid:20) cosh ( π py ) cosh ( π p γ / ) − pn tanh ( π p γ ) cosh ( π ny ) sinh ( π n γ / ) + p − n n (cid:21) . (C3)32he dispersion term < u z ϕ > in Eq. (43) can now be evaluated leading finally to the dispersionequation Eq. (31) with D eff given by Eq. (32). The numerical prefactor α is given by:1 α = − γπ ∞ ∑ n = , odd ∞ ∑ p = , even ∞ ∑ q = , even n p q ( p − n ) ( q − n ) (cid:20) n (cid:0) p tanh (cid:0) πγ p (cid:1) − q tanh (cid:0) πγ q (cid:1)(cid:1) p − q − np tanh (cid:0) πγ p (cid:1) (cid:0) n − q coth (cid:0) πγ n (cid:1) tanh (cid:0) πγ q (cid:1)(cid:1) n − q + ( p − n ) tanh (cid:0) πγ q (cid:1) q + ( p − n ) tanh (cid:0) πγ p (cid:1) p (cid:21) , (C4)and therefore only depends on the aspect ratio γ of the channel.Asymptotic approximations of α can be found for γ → α = γ + O ( γ ) . (C5)For a wide slit, i.e. γ (cid:29)
1, the terms tanh and coth in Eq. (C4) are close to (cid:39)
1, and one thus find:1 α (cid:39) − γ , (C6)where the factor (cid:39) α calculated using Eq. (C4) for several aspect ratios γ along with the asymp-totic behaviors Eqs. (C5) and (C6). It should be noted that the approximation Eq. (C6) valid for . -1 , . , = FIG. 15. Prefactor α in Eq. (32) for a rectangular cross-section of aspect ratio γ . The dashed dotted line isthe asymptotic bahavior Eq. (C5) valid for γ →
0. The magenta line is the approximation Eq. (C6) valid for γ (cid:29)
1. The dashed line is the value for a slit α = γ range 1–4. thin channels γ (cid:29) α within < .
5% for γ ≥ .
5. The value for asquare cross-section is α (cid:39) ppendix D: Derivation of the secondary transverse flow We introduce the stream function ψ ( x , y ) defined by: u x = ∂ ψ∂ y and u y = − ∂ ψ∂ x . (D1)Eqs. (52) and (53) turn to the inhomogeneous biharmonic equation: (cid:18) ∂ ∂ x + ∂∂ y (cid:19) ψ = Ra ∂ ϕ ∂ y , (D2)with ϕ given by Eq. (C3). The no-slip boundary conditions at the solid walls impose: ∂ ψ∂ y ( x , y = ± γ / ) = ∂ ψ∂ x ( x = , y ) = ψ ( x , y = ± γ / ) = ψ ( x = , y ) = . (D3)It can be seen that the same equations govern viscous flows induced by inhomogeneous tempera-ture fields in a rectangular container, and we will use the method described in Refs. [44, 45] to esti-mate the solution of Eq. (D2). In brief, the general solution of Eq. (D2) is written as ψ = ψ in + ψ in with ψ in any solution of the inhomogeneous equation (i.e. with the right-hand term), and ψ ho thesolution of the homogenous biharmonic equation, so that ψ fulfills the above boundary conditions. ψ in can be found simply by a Fourier series expansion following ϕ , and ψ ho can be found using adirect eigenfunction expansion, see below and Refs. [44, 45] for details.For simplicity, we re-write ϕ as: ϕ = Ra (cid:18) ∂ ϕ ∂ z (cid:19) ∞ ∑ n = , odd ∞ ∑ p = , even h n , p ( y ) cos ( n π x ) , (D4)with h n , p ( y ) given in Eq. (C3). An inhomogeneous solution ψ in of Eq. (D2) can be easily foundusing the following Fourier series: ψ in = (cid:18) Ra ∂ ϕ ∂ z (cid:19) ∞ ∑ n = , odd ∞ ∑ p = , even a n , p ( y ) cos ( n π x ) , (D5)with a n , p ( y ) solutions of the following ordinary differential equations: a ( ) n , p ( y ) − π n a ( ) n , p ( y ) + π n a n , p ( y ) = h ( ) n , p ( y ) . (D6)The solutions of these equations with the following boundary conditions: a ( ) n , p ( ± γ / ) = a n , p ( ± γ / ) = , (D7)34ead to a Fourier series representing ψ in . For the sake of clarity, the functions a n , p ( y ) are notwritten here. This Fourier series fullfills all the boundary conditions given by Eq. (D3) except ψ in ( x = , y ) = − ψ in ( x = , y ) = κ ( y ) (cid:54) = ψ ho of the homogeneous biharmonicequation with the boundary conditions given by Eq. (D3) but with ψ ho ( x = , y ) = − ψ ho ( x = , y ) = − κ ( y ) . As the symmetry of the problem imposes that ψ is an odd function of y , we willuse the following odd eigunfunctions exp ( ± λ n x ) φ n ( y ) with: φ n ( y ) = y cos ( λ n y ) − γ (cid:18) γλ n (cid:19) sin ( λ n y ) , (D8)and λ n the complex roots of the transcendental equation:sin ( λ n ) = λ n , (D9)which can be estimated by the Newton’s method. These eigenfunctions verify the boundary con-ditions expected at y = ± γ /
2. As ψ is a real and odd function of x with respect to x = /
2, aneigenfunction expansion for ψ ho ( x , y ) is: ψ ho ( x , y ) = (cid:18) Ra ∂ ϕ ∂ z (cid:19) ∞ ∑ n = R e (cid:104) b n φ n ( y ) (cid:16) e − λ n x − e − λ n ( − x ) (cid:17)(cid:105) (D10)where the complex numbers b n have to be determined from the boundary conditions at x = x =
1. To get an approximate solution, we proceed as proposed by Shankar in Ref. [44]. First, thesum in Eq. (D10) is truncated to the first N terms. Then, a least-squares procedure is used to findthe coefficients b n which yield the best expected boundary conditions at x = m equidistantpoints over the interval y = [ − γ / ] . This procedure corresponding to the resolution of 2 N linearalgebraic equations is performed using Mathematica. The coefficients b n are rapidly convergingand only a few eigenfunctions are needed to get an accurate estimate of ψ ho ( x , y ) .Figures 13(a) and 13(b) display the components u y and u x computed from the stream function ψ = ψ in + ψ ho calculated using the above procedure for γ =
1. Although the velocity field ( u x , u y )seem to suggest rotational symmetry for a square cross-section, this is not the case because of thegravity along x , and the maximum values of the components u x and u y , although close, are notstrictly equal. [1] N. Convery and N. Gadegaard, “30 years of microfluidics,” Micro and Nano Engineering , 76 (2019).
2] G. M. Whitesides, “The origins and the future of microfluidics,” Nature , 368 (2006).[3] D. J. Beebe, G. A. Mensing, and G. M. Walker, “Physics and applications of microfluidics in biology,”Annu. Rev. Biomed. Eng , 261 (2002).[4] H. A. Stone, A. D. Stroock, and A. Ajdari, “Engineering flows in small devices: Microfluidics Towarda Lab-on-a-Chip,” Annu. Rev. Fluid Mech. , 381 (2004).[5] T. M. Squires and S. R. Quake, “Microfluidics: fluid physics at the nanoliter scale,” Rev. Mod. Phys. , 977 (2005).[6] Y. Gu, V. Hegde, and K. J. M. Bishop, “Measurement and mitigation of free convection in microfluidicgradient generators,” Lab Chip , 3371 (2018).[7] S. K. Yoon, M. Mitchell, E. R. Choban, and P. J. A. Kenis, “Gravity-induced reorientation of theinterface between two liquids of different densities flowing laminarly through a microchannel,” LabChip , 1259 (2005).[8] B. Selva, L. Daubersies, and J.-B. Salmon, “Solutal convection in confined geometries: enhancementof colloidal transport,” Phys. Rev. Lett. , 198303 (2012).[9] A. M. J. Edwards, P. S. Atkinson, C. S. Cheung, H. Liang, D. J. Fairhurst, and F. F. Ouali, “Density-Driven Flows in Evaporating Binary Liquid Droplets,” Phys. Rev. Lett. , 184501 (2018).[10] Y. Li, C. Diddens, P. Lv, H. Wijshoff, M. Versluis, and D. Lohse, “Gravitational Effect in EvaporatingBinary Microdroplets,” Phys. Rev. Lett. , 114501 (2019).[11] K. H. Kang, H. C. Lim, H. W. Lee, and S. J. Lee, “Evaporation-induced saline Rayleigh convectioninside a colloidal droplet,” Phys. Fluids , 042001 (2013).[12] T. K. Pradhan and P. K. Panigrahi, “Convection inside condensing and evaporating droplets of aqueoussolution,” Soft Matter , 4335 (2018).[13] S. J. Lee, J. Hong, and Y.-S. Choi, “Evaporation-induced flows inside a confined droplet of dilutedsaline solution,” Langmuir , 7710 (2014).[14] L. Daubersies, J. Leng, and J.-B. Salmon, “Confined drying of a complex fluid drop: phase diagram,activity, and mutual diffusion coefficient,” Soft Matter , 5923 (2012).[15] C. Loussert, A. Bouchaudy, and J.-B. Salmon, “Drying dynamics of a charged colloidal dispersion ina confined drop,” Phys. Rev. Fluids , 084201 (2016).[16] K. Inoue and S. Inasawa, “Drying-induced back flow of colloidal suspensions confined in thin unidi-rectional drying cells,” RSC Advances , 15763 (2020).[17] T. K. Pradhan and P. K. Panigrahi, “Evaporation-induced natural convection of a liquid slug of binary ixture inside a microchannel: effect of confinement,” Microfluidics and Nanofluidics , 115 (2016).[18] R. Savino and R. Monti, “Buoyancy and surface-tension-driven convection in hanging-drop proteincrystallizer,” J. Crystal Growth , 308 (1996).[19] T. K. Pradhan, M. Asfer, and P. K. Panigrahi, “Droplet hydrodynamics during lysozyme proteincrystallization,” Phys. Rev. E , 051602 (2012).[20] V. Apostolopoulou, N. Junius, R. P. Sear, and M. Budayova-Spano, “Mixing Salts and Poly(ethyleneglycol) into Protein Solutions: The Effects of Diffusion across Semipermeable Membranes and ofConvection,” Cryst. Growth Des. , 3927 (2020).[21] T. K. Pradhan and P. K. Panigrahi, “Suppressing internal convection of a droplet using confinementduring protein crystallization,” J. Appl. Phys. , 084701 (2020).[22] J. Dunstan, K. J. Lee, Y. Hwang, S. F. Park, and R. E. Goldstein, “Evaporation-driven convectiveflows in suspensions of nonmotile bacteria,” Phys. Rev. Fluids , 123102 (2018).[23] C. L. Hansen, E. Skordalakes, J. M. Berger, and S. R. Quake, “A robust and scalable microfluidicmetering method that allows protein crystal growth by free interface diffusion,” Proc. Natl. Acad. Sci.USA , 16531 (2002).[24] A. Yamada, R. Renault, A. Chikina, B. Venzac, I. Pereiro, S. Coscoy, M. Verhulsel, M. C. Parrini,C. Villard, J.-L. Viovy, and S. Descroix, “Transient microfluidic compartmentalization using ac-tionable microfilaments for biochemical assays, cell culture and organs-on-chip,” Lab Chip , 4691(2016).[25] Y. Hallez and J. Magnaudet, “Effects of channel geometry on buoyancy-driven mixing,” Phys. Fluids , 053306 (2008).[26] J. O. Shin, S. B. Dalziel, and P. F. Linden, “Gravity currents produced by lock exchange,” J. FluidMech. , 1 (2004).[27] J. Martin, N. Rakotomalala, L. Talon, and D. Salin, “Viscous lock-exchange in rectangular channels,”J. Fluid Mech. , 132 (2011).[28] G. P. Matson and A. J. Hogg, “Viscous exchange flows,” Phys. Fluids , 023102 (2012).[29] T. S´eon, J. Znaien, D. Salin, J. P. Hulin E. J., Hinch, and B. Perrin, “Transient buoyancy-driven frontdynamics in nearly horizontal tubes,” Phys. Fluids , 123603 (2007).[30] M. L. Szulczewski and R. Juanes, “The evolution of miscible gravity currents in horizontal porouslayers,” J. Fluid Mech. , 82 (2013).[31] M. T. Tyn and W. F. Calus, “Temperature and Concentration Dependence of Mutual Diffusion Coeffi- ients of Some Binary Liquid Systems,” J. Chem. Eng. Data , 310 (1975).[32] J. Crank, The mathematics of diffusion (Oxford university press, 1975).[33] S. Xin and P. Le Qu´er´e, “An extended Chebyshev pseudo-spectral benchmark for the 8:1 differentiallyheated cavity,” Int. J. Heat Mass Transfer , 981 (2002).[34] S. Xin, J. Chergui, and P. Le Qu´er´e, “3D spectral parallel multi-domain computing for natural con-vection flows,” in Parallel Computational Fluid Dynamics (2008) pp. 163–171.[35] M. E. Erdogan and P. C. Chatwin, “The effects of curvature and buoyancy on the laminar dispersionof solute in a horizontal tube,” J. Fluid Mech. , 465 (1967).[36] R. Smith, “Longitudinal dispersion of a buoyant contaminant in a shallow channel,” J. Fluid Mech. , 677 (1976).[37] N. G. Barton, “The dispersion of a buoyant solute in laminar flow in a straight horizontal pipe. Part 1.Predictions from Erdogan & Chatwin’s (1967) paper,” J. Fluid Mech. , 81 (1976).[38] J.S. Godfrey, “A numerical model of the James River estuary, Virginia, U.S.A.” Estuarine Coastal Mar.Sci. , 295 (1980).[39] D. J. Maclean and T. Alboussi`ere, “Measurement of solute diffusivities. Part I. Analysis of coupledsolute buoyancy-driven convection and mass transport,” Int. J. Heat Mass Transfer , 1639 (2001).[40] J.-B. Salmon and F. Doumenc, “Buoyancy-driven dispersion in confined drying of liquid binary mix-tures,” Phys. Rev. Fluids , 024201 (2020).[41] W. R. Young and Scott Jones, “Shear dispersion,” Phys. Fluids A , 1087 (1991).[42] P. C. Chatwin and P. J. Sullivan, “The effect of aspect ratio on longitudinal diffusivity in rectangularchannels,” J. Fluid Mech. , 347 (1982).[43] A. Ajdari, N. Bontoux, and H. A. Stone, “Hydrodynamic dispersion in shallow microchannels: theeffect of cross-sectional shape,” Anal. Chem. , 387 (2006).[44] P. N. Shankar, “The eddy structure in Stokes flow in a cavity,” J. Fluid Mech. , 371 (1993).[45] P. N. Shankar, V. V. Meleshko, and E. I. Nikiforovich, “Slow mixed convection in rectangular con-tainers,” J. Fluid Mech. , 203 (2002).[46] D. D. Joseph, A. Huang, and H. Hu, “Non-solenoidal velocity effects and Korteweg stresses in simplemixtures of incompressible liquids,” Physica D , 104 (1996).[47] H. Brenner, “Navier-Stokes revisited,” Physica A , 60 (2005).[48] H. Bruus, Theoretical Microfluidics , edited by Oxford Master Series in Physics (2007)., edited by Oxford Master Series in Physics (2007).