# Role of solutal free convection on interdiffusion
in a horizontal microfluidic channel

Jean-Baptiste Salmon  CNRS, Solvay, LOF, UMR 5258, Univ. Bordeaux, F-33600
Pessac, France.  Laurent Soucasse  Laboratoire EM2C, CNRS, CentraleSupélec,
Université Paris-Saclay, Grande Voie des Vignes, 92295 Chatenay-Malabry cedex,
France.  Frédéric Doumenc  Université Paris-Saclay, CNRS, FAST, 91405, Orsay,
France,  Sorbonne Université, UFR 919, 4 place Jussieu, F-75252, Paris Cedex
05, France.

###### Abstract

We theoretically investigate the role of solutal free convection on the
diffusion of a buoyant solute at the microfluidic scales, $\simeq
5$–$500~{}\mu$m. We first consider a horizontal microfluidic slit, one half of
which initially filled with a binary solution (solute and solvent), and the
other half with pure solvent. The buoyant forces generate a gravity current
that couples to the diffusion of the solute. We perform numerical resolutions
of the 2D model describing the transport of the solute in the slit. This study
allows us to highlight different regimes as a function of a single parameter,
the Rayleigh number Ra which compares gravity-induced advection to solute
diffusion. We then derive asymptotic analytical solutions to quantify the
width of the mixing zone as a function of time in each regime and establish a
diagram that makes it possible to identify the range of Ra and times for which
buoyancy does not impact diffusion. In a second step, we present numerical
resolutions of the same model but for a 3D microfluidic channel with a square
cross-section. We observe the same regimes as in the 2D case, and focus on the
dispersion regime at long time scales. We then derive the expression of the 1D
dispersion coefficient for a channel with a rectangular section, and analyse
the role of the transverse flow in the particular case of a square section.
Finally, we show that the impact of this transverse flow on the solute
transport can be neglected for most of the microfluidic experimental
configurations.

###### pacs:

##  I Introduction

Microfluidics refers to a wide range of technological tools for manipulating
liquids in microfabricated networks of channels with cross-sectional
dimensions ranging from a few microns to a few hundred microns. Applications
of this technology are numerous and diverse, from high throughput miniaturized
bioassays to fundamental studies in physical-chemistry, see Refs. Convery and
Gadegaard (2019); Whitesides (2006); Beebe  _et al._ (2002) for some reviews.

Figure 1: (a) Schematic diagram of the diffusion experiment in a microfluidic
channel of height $H$ and width $L$. The red arrows represent the gravity
current $U_{\mathrm{B}}$ induced by the difference in density. The color code
indicates the concentration of the solute. (b) Width of the mixing zone $W$
vs. time $T$. The dashed line is the case for which buoyancy does not affect
the spreading of the solute.

The small scales of microfluidic technologies allow to study numerous
processes while controlling finely all transport phenomena (mass, momentum,
energy) Stone  _et al._ (2004). In particular, the role of buoyancy has been
mentioned by Squires and Quake’s in their review on the physics of fluids at
the nanoliter scale and quantified using scaling arguments as explained below
Squires and Quake (2005). 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 $\Phi_{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\sim\sqrt{DT}$ where $D$ is the diffusion coefficient
of the mixture and $T$ the time. Now assuming that the density evolves
linearly with the concentration in solute, i.e.:

$\displaystyle\rho=\rho_{0}(1+\beta\Phi)\,,$ (1)

where $\rho_{0}$ is the density of the solvent, buoyancy induces a
gravitational current, the solution flowing under the solvent for $\beta>0$,
see Fig. 1(a). The magnitude of the gravity current $U_{\mathrm{B}}$ can be
estimated from a balance between buoyant forces $\sim\rho_{0}\beta\Phi_{i}g$
and viscous forces $\sim\rho_{0}\nu U_{\mathrm{B}}/H^{2}$, leading to the
velocity scale:

$\displaystyle U_{\mathrm{B}}\sim\frac{\beta\Phi_{i}gH^{2}}{\nu}\,,$ (2)

where $\nu$ is the kinematic viscosity of the mixture Squires and Quake
(2005). Note that such a flow exists whatever the height $H$ of the channel as
the density gradient is orthogonal to the gravity field $\mathbf{g}$. The
impact of this flow on the solute transport can be determined using the
Rayleigh number:

$\displaystyle\text{Ra}=\frac{\beta\Phi_{i}gH^{3}}{\nu
D}\sim\frac{U_{\mathrm{B}}H}{D}\,,$ (3)

comparing advection to diffusion. $\mathrm{Ra}$ is similar to the Péclet
number describing mass transport in forced convection problems, but with the
important difference that the solute in the present case is not passive, since
the solute itself is the source of the flow Squires and Quake (2005).
Furthermore, a balance between viscous forces $\sim\rho_{0}\nu
U_{\mathrm{B}}/H^{2}$ and inertial forces $\sim\rho_{0}U_{\mathrm{B}}^{2}/H$
leads to the definition of the Grashof number, similarly to the Reynolds
number for forced convection:

$\displaystyle\text{Gr}=\frac{\beta\Phi_{i}gH^{3}}{\nu^{2}}=\frac{{\rm
Ra}}{{\rm Sc}}\,,$ (4)

where $\mathrm{Sc}=\nu/D$ is the Schmidt number. As ${\rm Sc}\geq 10^{2}$ for
most liquid mixtures, one has thus $\mathrm{Gr}\ll{\rm Ra}$ and viscous
dissipation a priori dominates inertia within the gravity current in most
microfluidic applications Squires and Quake (2005).

Many groups have reported gravity-driven currents in microfluidic experiments
for which density gradients are imposed either using membranes Gu  _et al._
(2018) or transverse mixing between coflowing miscible liquids Yoon  _et al._
(2005); Selva  _et al._ (2012). In a different context, many groups also
reported such flows when density gradients are induced by the evaporation of a
liquid mixture in a confined geometry ($H=100$–$1000~{}\mu$m) such as sessile
drops Edwards  _et al._ (2018); Li  _et al._ (2019); Kang  _et al._ (2013),
confined drops Pradhan and Panigrahi (2018); Lee  _et al._ (2014); Daubersies
_et al._ (2012); Loussert  _et al._ (2016), or micro-capillaries Inoue and
Inasawa (2020); Pradhan and Panigrahi (2016). The impact of buoyancy on the
solute transport is not always mentioned in such works, and most groups
consider that solutal free convection plays little role at the microfluidic
scales although experimental configurations with density gradients are
ubiquitous in applications.

To illustrate this point, let us consider interdiffusion between an aqueous
NaCl solution at a concentration of 1 M and pure water. For such a mixture,
$D\simeq 1.6\times 10^{-9}$ m2/s, the difference in density is $\simeq 38$
kg/m3 leading to $\beta\Phi_{i}\simeq 3.8\times 10^{-2}$, and the kinematic
viscosity is $\nu\simeq 10^{-6}~{}$m2/s Pradhan and Panigrahi (2016). For a
microfluidic channel of height $H=20~{}\mu$m, one finds $\text{Ra}\simeq 1.9$,
but due to the scaling $\mathrm{Ra}\propto H^{3}$, the Rayleigh number
increases to $\text{Ra}\simeq 230$ for $H=100~{}\mu$m, and even up to
$\text{Ra}\simeq 2.9\times 10^{4}$ for $H=500~{}\mu$m. This numerical
application illustrates the importance of going beyond the scaling laws
presented above to quantitatively predict the range of Rayleigh numbers for
which free convection plays only a minor role in a microfluidic configuration.
Furthermore, although these buoyancy-induced flows could have little impact on
the solute concentration gradients that generate them, they still do exist,
and are able to effectively disperse less mobile species in the fluid mixture,
such as macromolecules or colloids Selva  _et al._ (2012). These buoyancy-
driven flows may also have an influence in protein crystallization experiments
Savino and Monti (1996); Pradhan  _et al._ (2012); Apostolopoulou  _et al._
(2020); Pradhan and Panigrahi (2020), for evaluating colloidal diffusio-
phoresis induced by solute gradients Gu  _et al._ (2018), or even in the
context of biological systems for the motility of microorganisms Dunstan  _et
al._ (2018). It is therefore necessary to quantify these flows as a function
of the density gradients that generate them for predicting their possible role
on other species in the case of complex fluid mixtures.

In the present work, we study in depth the configuration presented in Fig. 1
in a microfluidic context, i.e. height $H$ in the $\simeq 5$–$500~{}\mu$m
range. This experimental configuration has been implemented many times in
microfluidic devices, either using valves or sliding walls for instance, for
various applications such as protein crystallization or biochemical assays
Hansen  _et al._ (2002); Yamada  _et al._ (2016). Our main goal is to
quantitatively delineate the range of Rayleigh numbers for which mixing is
impacted by buoyancy in such a configuration, and to predict the laws $W$ vs.
$T$. Configurations similar to that shown in Fig. 1 have been studied for
chemical or civil engineering applications and environmental issues that
involve length scales $H$ ranging typically from $0.1$ to $100$ m. In such
cases, commonly referred to in the literature as the ”lock-exchange” problem,
molecular diffusion is negligible and inertial effects are often significant
Hallez and Magnaudet (2008); Shin  _et al._ (2004); Martin  _et al._ (2011);
Matson and Hogg (2012); Séon  _et al._ (2007). Special mention should be made
however of the work of Szulczewski and Juanes Szulczewski and Juanes (2013),
who studied a situation very similar to that shown in Fig. 1(a) including also
molecular diffusion, but for a 2D porous layer in the context of geological
CO2 sequestration. Surprisingly, we are not aware of any work that has studied
the microfluidic case where diffusion cannot be neglected, which also
motivated this work. Since our work is related to microfluidic applications,
we have explored Rayleigh numbers up to $\mathrm{Ra}=10^{5}$. These high
values are at the limit of most microfluidic dimensions, but can easily be
obtained as soon as $H$ reaches the millimeter scale, even for dilute
solutions. For example, in the numerical application given previously, water
and salty water at 1 M, $\mathrm{Ra}=10^{5}$ for $H\simeq 750~{}\mu$m. As we
subsequently consider liquid mixtures only, the smallest Schmidt number we
explored is $\mathrm{Sc}=10^{2}$. Such small values can be observed in the
case of the diffusion of small molecules in a low-viscosity solvent, e.g.
water in acetone Tyn and Calus (1975).

The present paper is organized as follows. In Sec. II, we present the set of
equations modeling the transport of the solute in the configuration shown in
Fig. 1, as well as details about the numerical resolutions. In Sec. III, we
study the case of a 2D slit for the sake of simplicity, i.e. two infinite and
parallel plates separated by a thickness $H$. The numerical data show a rich
temporal succession of different regimes of solute spreading, that can be
captured using analytical asymptotic solutions, and the analogy with the case
of a 2D porous layer Szulczewski and Juanes (2013) is discussed. We finally
address in Sec. IV the case of a 3D microfluidic channel for which transverse
flows also exist. We finally conclude our work in Sec. V and insist on its
possible implications.

##  II Model and numerical resolution

###  II.1 Model and dimensionless variables

We consider the situation described in Fig. 1(a): a straight microfluidic
channel of infinite length and rectangular cross-section initially filled with
a solution at concentration $\Phi_{i}$ for $Z>0$ and by pure solvent for
$Z<0$. For the sake of simplicity, we consider that the kinematic viscosity
$\nu$ and the interdiffusion coefficient $D$ are constant, and that the
density of the solution evolves linearly with the volume fraction in solute,
Eq. (1). Because of this linearity, the problem described here 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 governing the solute transport and the velocity field $\mathbf{U}$
are:

$\displaystyle\rho_{0}\left(\frac{\partial\mathbf{U}}{\partial
T}+\mathbf{U}.\nabla\mathbf{U}\right)=\rho_{0}\nu\Delta\mathbf{U}-\nabla
P+(\rho(\Phi)-\rho_{0})\mathbf{g},$ (5) $\displaystyle\nabla.\mathbf{U}=0,$
(6) $\displaystyle\frac{\partial\Phi}{\partial
T}+\mathbf{U}.\nabla\Phi=D\Delta\Phi,$ (7)

where $P$ is the pressure deviation from the hydrostatic pressure field for
$\Phi=0$. Boundary conditions at the solid walls are the no-slip and the
impermeability conditions, $\mathbf{U}=0$ for the velocity field and
$\mathbf{n}.\nabla\Phi=0$ for the concentration field. We also impose
$\mathbf{U}(Z\to\pm\infty)=0$ and $\partial_{Z}\Phi(Z\to\pm\infty)=0$
resulting in no pressure-driven flow along the channel (i.e. the gravity
current is the only flow). Initial conditions are given by $\mathbf{U}=0$,
$\Phi=\Phi_{i}$ for $Z>0$ and $\Phi=0$ for $Z<0$.

To get more insights into the mechanisms of solute transport, we define the
following dimensionless variables:

$\displaystyle
x=X/H,~{}~{}y=Y/H,~{}~{}z=Z/H,~{}~{}t=DT/H^{2},~{}~{}\gamma=L/H\,,$ (8)
$\displaystyle\mathbf{u}=(H/D)\mathbf{U},~{}~{}p=H^{2}/(\rho_{0}\nu
D)P,~{}~{}\varphi=\Phi/\Phi_{i}\,.$ (9)

With such definitions, the model given by Eqs. (5-7) reads now:

$\displaystyle\frac{1}{\text{Sc}}\left(\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}.\nabla\mathbf{u}\right)=\Delta\mathbf{u}-\nabla
p-\text{Ra}\varphi\mathbf{e_{x}}\,,$ (10)
$\displaystyle\nabla.\mathbf{u}=0\,,$ (11)
$\displaystyle\frac{\partial\varphi}{\partial
t}+\mathbf{u}.\nabla\varphi=\Delta\varphi\,,$ (12)

where $\mathbf{e_{x}}$ is the unit vector along $x$. With these variables, the
initial conditions are:

$\displaystyle\mathbf{u}(x,y,z,t=0)=\mathbf{0}~{}~{}{\rm
and}~{}~{}\varphi(x,y,z,t=0)=\mathcal{H}(z)\,,$ (13)

where $\mathcal{H}(z)$ is the Heaviside function. Boundary conditions are
given by:

$\displaystyle\mathbf{u}=0~{}~{}{\rm and}~{}~{}\mathbf{n}.\nabla\varphi=0$
(14)

on the solid walls, and:

$\displaystyle\mathbf{u}(z\to\pm\infty)=0~{}~{}{\rm
and}~{}~{}\partial_{z}\varphi(z\to\pm\infty)=0\,.$ (15)

To estimate the role of buoyancy on the solute spreading, we first define the
cross-section averaged concentration profile by:

$\displaystyle\varphi_{0}(z,t)=<\varphi>=\frac{1}{\gamma}\int_{-\gamma/2}^{\gamma/2}\int_{0}^{1}\varphi(x,y,z,t)\text{d}x\text{d}y\,,$
(16)

and the extent of the mixing zone by:

$\displaystyle
w(t)=\sqrt{\frac{1}{2}\int_{-\infty}^{\infty}{z}^{2}\frac{\partial\varphi_{0}}{\partial
z}\text{d}z}\,.$ (17)

In the case of a neutrally-buoyant solute, $\mathrm{Ra}=0$ and the model
described above admits the following simple solution Crank (1975):

$\displaystyle\mathbf{u}=0\,,$ (18)
$\displaystyle\varphi(x,y,z,t)=\varphi_{0}(z,t)=\frac{1}{2}\left[1+\text{Erf}\left(\frac{z}{2\sqrt{t}}\right)\right]\,.$
(19)

In this case, the width of the interdiffusion zone is given by $w=\sqrt{t}$,
i.e. $W=\sqrt{DT}$ with real units. This is the classical square-root
spreading of the solute due to molecular diffusion. For $\mathrm{Ra}>0$, any
deviation from this simple law is a priori a signature of buoyancy-induced
dispersion, see Fig. 1(b).

###  II.2 Numerical resolution

Eqs. (10-14) have been solved numerically for two distinct geometries: the 2D
case of a microfluidic slit ($\gamma\to\infty$) and the 3D case of a
rectangular micro-channel with a square cross-section ($\gamma=1$). In both
cases, the boundary conditions Eq. (15) at $z\to\pm\infty$ have been moved to
$z=\pm\lambda$ where $\lambda$ is a finite distance such that $\lambda\gg w$.

The 2D numerical simulations have been performed with the commercial software
Comsol Multiphysics based on finite elements (Galerkin method). Time
discretization is based on implicit Backward Differentiation Formulas, with an
adaptive time stepping. Spatial discretization was achieved by a structured
mesh of Lagrangian elements, linear for the pressure and quadratic for the
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$ GB of RAM.

The 3D simulations required the use of in-house made software Xin and Le Quéré
(2002); Xin  _et al._ (2008), specifically optimized for the simulation of
free convection in cavities on parallel architectures and based on a
multidomain spectral method. Chebyshev collocation is used for spatial
discretization of the three dimensions of space. The pressure-flow coupling is
ensured by a projection method that forces the velocity divergence-free
condition. Time integration is performed through a second order temporal
scheme combining a Backward Differentiation (BDF2) scheme for the linear terms
with an Adams Bashforth extrapolation for the convective terms. Domain
decomposition along the $z$-horizontal direction is carried out by the Schur
complement method for parallelization purposes. Each spatial domain is a cube
of size $1$ in dimensionless units. The spatial resolution in the
$z$-direction has been increased at the first moments of the simulation in
order to capture the stiff concentration gradient and the mesh convergence has
been checked by observing the decay of the Chebyshev spectral coefficients.
Computations were made in a HPC facility, using from 40 to 360 Intel Xeon 2.30
GHz processors.

For both 2D and 3D geometries, simulations have been divided into several time
intervals in order to adapt the simulation parameters to the temporal
evolution of the flow. For each new time interval, the length $2\lambda$ of
the spatial domain was extended to take into account the increase of the
mixing zone.

##  III 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 fixed Schmidt number $\text{Sc}=10^{5}$,
and three different Rayleigh numbers $\text{Ra}=10^{3}$, $10^{4}$, and
$10^{5}$, over a wide range of time scales, from $t=10^{-7}$ to $t\simeq
10^{6}$. To test the role of the Schmidt number, we also performed numerical
resolutions for the same Rayleigh numbers, $\mathrm{Sc}=10^{2}$, $10^{3}$,
$10^{4}$ and $10^{5}$, and time scales ranging from $t=10^{-7}$ to $t\simeq
5\times 10^{-3}$ as $\mathrm{Sc}$ only plays a role at early time scales, see
below.

###  III.1 The case $\text{Ra}=10^{5}$ and $\text{Sc}=10^{5}$

We begin with the case $\text{Ra}=10^{5}$ and $\text{Sc}=10^{5}$. Figure 2(a)
displays several 2D concentration maps $\varphi(x,z,t)$ at the times shown in
Fig. 2(b), see also movie M1 corresponding to these data in the ESI.

Figure 2: For $\text{Ra}=10^{5}$ and $\text{Sc}=10^{5}$: (a) concentration
fields $\varphi(x,z,t)$ at the times indicated by the corresponding symbols in
(b), the $z$-scale is different in each image. (b) Width of the mixing zone
$w(t)$. The magenta dashed line is given by Eq. (27) and corresponds to early
advection with $w\sim t$. The yellow dashed line is given by Eq. (29) and
corresponds to late advection with $w\sim\sqrt{t}$. The green dashed line is
computed from the numerical resolution of Eq. (31) and includes both the 1D
dispersion regime $w\sim t^{1/4}$ and late diffusion $w\sim\sqrt{t}$. The blue
dashed line is the diffusion law $w=\sqrt{t}$. See also movie M1 corresponding
to these data in the ESI.

These data evidence a combination of solute spreading by diffusion and
buoyancy-induced advection.

Figure 2(b) displays $w(t)$ computed from the numerical resolution for the
case $\mathrm{Ra}=10^{5}$ and $\mathrm{Sc}=10^{5}$, along with the diffusion
law $w=\sqrt{t}$. These data clearly show that the mixing zone evolves
according to the diffusion law expected without buoyancy at small time scales
$t\ll 10^{-5}$, but also at long time scales $t\gg 10^{4}$. For intermediate
times, the width of the mixing zone is significantly larger, evidencing the
role of the buoyancy-driven advection. The $w$ vs. $t$ behavior can be
rationalized using different regimes, each with a given power law $w\sim
t^{\delta}$ shown in Fig. 2(b). These regimes are presented below in detail,
along with self-similar asymptotic solutions for the concentration profiles
$\varphi_{0}(z,t)$ and the corresponding spreading laws $w$ vs. $t$.

###  III.2 Early diffusion regime

We first analyse the transport of solute at short time scales.

Figure 3: (a) Height-averaged concentration profiles $\varphi_{0}(z,t)$ vs.
$z$ for time scales ranging from $t=10^{-7}$ to $t\leq 1.6\times 10^{-6}$ (7
curves, $\text{Ra}=10^{5}$ and $\text{Sc}=10^{5}$). (b) Same data plotted
against the reduced variable $z/\sqrt{t}$, the dashed line is given by Eq.
(19).

Figure 3(a) displays the height-averaged concentration profile
$\varphi_{0}(z,t)$ [Eq. (16)] for time scales $t\leq 1.6\times 10^{-6}$
evidencing the spreading of the solute. As shown in Fig. 3(b), all the
profiles collapse on a single curve when plotted against the reduced variable
$z/\sqrt{t}$. This curve is correctly described by Eq. (19) demonstrating 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 computed from the 2D data using Eq. (17)
is correctly fitted by $w=\sqrt{t}$ for $t\leq 1.6\times 10^{-6}$ as shown in
Fig. 2(b).

However, a flow exists as $\mathbf{u}=0$ does not satisfy the Navier-Stokes
equation Eq. (10), because of the non zero $z$-component of the pressure
gradient due to buoyancy. Figure 4(a) indeed displays the flow driven by this
difference of density for $t\simeq 1.3\times 10^{-6}$. This gravity current
corresponds to a recirculating flow developed within the slit on a length
scale $\sim 1$. Figure 4(b) showing the maximal value of the component $u_{z}$
in the plane $z=0$ in this early regime, evidences that the velocity steadily
increases up to reaching a plateau value of $\simeq 10^{3}$ for time scales
$t\geq 0.7\times 10^{-6}$ for $\mathrm{Sc}=10^{5}$. This plot also shows the
same data corresponding to several Schmidt numbers, $\mathrm{Sc}=10^{2}$,
$10^{3}$, $10^{4}$ and $10^{5}$. After a transient, the velocity reaches the
same plateau value for all Schmidt numbers, except for $\mathrm{Sc}=10^{2}$
for which the plateau is not reached.

Figure 4: (a) Velocity vector field superimposed with the concentration field
at $t=1.3\times 10^{-6}$ ($\mathrm{Ra}=10^{5}$, $\mathrm{Sc}=10^{5}$). (b)
Maximal value of the component $u_{z}$ in the plane $z=0$ vs. $t$ for
$\mathrm{Ra}=10^{5}$ and $\mathrm{Sc}=10^{2}$ ($\star$), $10^{3}$ ($\square$),
$10^{4}$ ($\circ$), $10^{5}$ ($\triangledown$). The horizontal dashed line is
$\simeq 1060$. (c) Maximal value of $u_{z}(x,z=0)$ rescaled by $\mathrm{Ra}$
vs. $\mathrm{Sc}\,t$ for $\mathrm{Ra}=10^{5}$ (blue), $10^{4}$ (magenta), and
$10^{3}$ (cyan) and $\mathrm{Sc}=10^{2}$ ($\star$), $10^{3}$ ($\square$),
$10^{4}$ ($\circ$), $10^{5}$ ($\triangledown$). The horizontal dashed line is
estimated using Eq. (22) and given by $\simeq 0.0106$, the vertical dashed
line is $\mathrm{Sc}\,t=0.0571$.

As the transients depend on $\mathrm{Sc}$ and thus on the Grashof number
$\mathrm{Gr}$, see Eq. (4), these results suggest the existence of an inertial
regime corresponding to the development of the recirculating gravity current
through the slit. To better highlight this regime, Fig. 4(c) displays the
maximal $z$-component of the velocity field at $z=0$ vs. time $t$ for several
Rayleigh numbers $\mathrm{Ra}=10^{3}$, $10^{4}$, $10^{5}$ and several Schmidt
numbers $\mathrm{Sc}=10^{2}$, $10^{3}$, $10^{4}$ and $10^{5}$. All the
transients collapse on a single curve when times are scaled by $1/\text{Sc}$
and velocities by $\mathrm{Ra}$. This result is recovered from the Navier-
Stokes equation Eq. (10) assuming that the non-linear inertial term
$\mathbf{u}.\nabla\mathbf{u}$ does not play any role. The start-up of the flow
therefore corresponds simply to the diffusion of the momentum through the
slit, expected to take place on a time scaling as $\sim 1/\text{Sc}$. By
estimating the time it takes for the maximal velocity $u_{z}$ to reach 90% of
its plateau value, the duration of the inertial regime is $t\simeq
0.0571/\text{Sc}$, see the vertical dashed line in Fig. 4(c). The smallest
Schmidt number investigated, $\mathrm{Sc}=10^{2}$, is an exception. In this
case, the steady plateau is not observed because inertia is still significant
after the end of the diffusion regime, i.e. when advection starts affecting
the solute transport. Turning to real units, the duration of the inertial
regime is given by $\simeq 0.0571\,H^{2}/\nu$, and lasts only a few tens of
milliseconds even for $H=500~{}\mu$m and low-viscosity solvents $\nu=5\times
10^{-7}$ m2/s. This numerical application shows that such a regime cannot be
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 transport of the solute.

We now turn to the part of the diffusive regime characterized by a steady
gravity current, i.e. 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, and
it remains close to the initial condition Eq. (13), as evidenced by Fig. 4(a)
for the case $\mathrm{Ra}=10^{5}$. The velocity field after the inertial
transient is therefore expected to be the solution of the steady Stokes
equation:

$\displaystyle 0=\Delta\mathbf{u}-\nabla
p-\text{Ra}\mathcal{H}(z)\mathbf{e_{x}}\,,$ (20)

where $\nabla.\mathbf{u}=0$ and $\mathcal{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=1/2$ and
$u_{z}$ in the plane $z=0$:

$\displaystyle
u^{p}_{x}(x=1/2,z)=\frac{2\mathrm{Ra}}{\pi}\int_{0}^{\infty}\frac{\sinh^{2}\left(\frac{k}{4}\right)\left[k-2\sinh\left(\frac{k}{2}\right)\right]\sin(kz)}{k^{3}[k+\sinh(k)]}\,\text{d}k\,,$
(21) $\displaystyle
u^{p}_{z}(x,z=0)=-\frac{\mathrm{Ra}}{\pi}\int_{0}^{\infty}\frac{(x-1)\sinh(kx)+x\sinh(k-kx)}{k^{2}[k+\sinh(k)]}\,\text{d}k\,.$
(22)

These two expressions correctly approximate the velocity profiles in the
plateau regime, see the black lines in Figs. 5(a) and 5(b) for the case
$\mathrm{Ra}=10^{5}$ and $\mathrm{Sc}=10^{5}$. The maximal velocities are
about $\text{max}[u^{p}_{z}(x,z=0)]\simeq 0.0106\,\text{Ra}$ [at $x\simeq
0.81$, the plateau value in Fig. 4(c)], and
$\text{max}[u^{p}_{x}(x=1/2,z)]\simeq 0.0065\,\text{Ra}$ (at $z\simeq-0.23$).

Figure 5: (a) $u_{x}(x=1/2,z)$ vs. $z$ and (b) $u_{z}(x,z=0)$ vs. $x$ for time
scales ranging from $t=1.6\times 10^{-6}$ to $10^{-4}$ (10 curves,
$\mathrm{Ra}=10^{5}$ and $\mathrm{Sc}=10^{5}$). The thin black lines are given
by Eq. (21) in (a) and by Eq. (22) in (b).

###  III.3 Early advection

The initial diffusion regime described in Sec. III.2 ceases when the effect of
advection on the the solute transport becomes non negligible. This transition
time corresponds to the departure from the diffusion law $w=\sqrt{t}$ at
$t\simeq 1\times 10^{-5}$ for the case $\mathrm{Ra}=10^{5}$ and
$\mathrm{Sc}=10^{5}$, see Fig. 2(b). A typical 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:

Figure 6: Concentration fields $\varphi$ at times (a) $t\simeq 6.3\times
10^{-5}$ and (b) $t\simeq 6.3\times 10^{-4}$ for $\mathrm{Ra}=10^{5}$ and
$\mathrm{Sc}=10^{5}$. The magenta dashed lines are the deformations estimated
by Eq. (25). (c) Average concentration profiles against $z/(\mathrm{Ra}\,t)$
for time scales ranging from $t\simeq 2\times 10^{-5}$ to $1.3\times 10^{-4}$
(9 curves), the dashed line is given by Eq. (26). Inset: same profiles against
$z$.

####

The 2D deformation of the concentration field, obviously due to advection, is
much larger than the diffusive spreading. For this reason, we can neglect
diffusion against advection and approximate the concentration field by an
Heaviside function:

$\varphi(x,z,t)=\mathcal{H}[z-z_{f}(x,t)]\,,$ (23)

where $z_{f}(x,t)$ is the position of a front separating two regions, one
where $\varphi\simeq 0$ and the other one where $\varphi\simeq 1$. In this
advection regime, the front motion along $z$ during a small time interval
$\mathrm{d}t$ reads:

$\mathrm{d}z_{f}(x,t)=u_{z}(x,z_{f}(t),t)\,\mathrm{d}t\,.$ (24)

####

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 given by Eqs. (21) and (22) corresponding to
the plateau observed in Fig. 4(b). We deduce from these assumptions and Eq.
(24) that the front profile can be approximated by:

$z_{f}(x,t)\simeq u_{z}^{p}(x,z=0)\,t\,,$ (25)

where $u_{z}^{p}(x,z=0)$ is the steady velocity field given by Eq. (22). Eq.
(25) is consistent with the concentration 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 in this regime:

$\varphi_{0}(z,t)=\int_{0}^{1}\mathcal{H}[z-u_{z}^{p}(x,z=0)\,t]\,\mathrm{d}x\,.$
(26)

One can easily show that the above relation is a self-similar function of the
variable $z/(\mathrm{Ra}\,t)$. Figure 6(c) compares the theoretical relation
Eq. (26) with the numerical simulations, evidencing a reasonably good collapse
of the data on the theoretical master curve. Furthermore, the width of the
mixing zone $w(t)$ defined by Eq. (17) can be estimated using Eq. (26),
leading to:

$w\simeq 0.00538\mathrm{Ra}\,t\,.$ (27)

This behavior is plotted in Fig. 2(b) and correctly fits the data obtained
from the numerical resolution of the 2D model from $t=1.3\times 10^{-5}$ to
$2\times 10^{-4}$.

At 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=6.3\times 10^{-4}$. Indeed, the deformation of the
concentration field is of the order of $1$ for such time scales and Eq. (22)
can no longer be used for the estimation of the velocity field. It marks the
end of the early advection regime.

###  III.4 Late advection

At later time scales, solute spreading by diffusion still remains negligible
as compared to solute advection, but the 2D deformation of the concentration
fields is now much larger than $1$, see for instance Fig. 7(a) showing a
snapshot at $t=1.6\times 10^{-2}$.

Figure 7: (a) Concentration field $\varphi$, and components (b) $u_{z}$, (c)
$u_{x}$ at $t=1.6\times 10^{-2}$ for $\mathrm{Ra}=10^{5}$ and
$\mathrm{Sc}=10^{5}$. The magenta dashed line is given by the solution of Eq.
(28) in the case of a slit. (c) Average concentration profiles
$\varphi_{0}(z,t)$ vs. $\xi=z/\sqrt{\tilde{D}t}$ for time scales ranging from
$t\simeq 0.003$ to $0.02$ (9 curves), the dashed line is $\psi(\xi)$ solution
of Eq. (28). Inset: same profiles against $z$.

These data along with the flow field shown in Figs. 7(b) and 7(c), evidence
the reciprocal exchange of the solution and the solvent separated by a diffuse
pseudo-interface. For fully negligible diffusion, this regime commonly
referred to as the viscous 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\sim\sqrt{\tilde{D}t}$ where
$\tilde{D}$ is an effective diffusion coefficient. The square-root behavior
arises from the competition between buoyant forces
($\sim\rho_{0}\beta\Phi_{i}gH/W$) and viscous forces
($\sim\rho_{0}\nu\dot{W}/H^{2}$), leading to the dimensionless scaling law
$w^{2}\sim\mathrm{Ra}\,t$ Matson and Hogg (2012). The effective diffusion
coefficient $\tilde{D}$ therefore only depends on $\mathrm{Ra}$ and the
geometry, and it has been calculated in various cases: porous medium, circular
tube, but also rectangular channel and slit Szulczewski and Juanes (2013);
Séon  _et al._ (2007); Matson and Hogg (2012); Martin  _et al._ (2011). The
main idea of these works is to compute the shape $\psi(z,t)$ of the pseudo-
interface separating the two fluids, assuming large deformations and thus a
quasi-parallel flow along the $z$-axis (lubrication approximation). With such
approximations, one can show that $\psi(z,t)$ admits a self-similar shape
$\psi(\xi)$ with $\xi=z/\sqrt{\tilde{D}t}$, solution of:

$\displaystyle-\xi\frac{\mathrm{d}\psi}{\mathrm{d}\xi}=2\frac{\mathrm{d}}{\mathrm{d}\xi}\left(f(\psi)\frac{\mathrm{d}\psi}{\mathrm{d}\xi}\right)\,,$
(28)

with $f(\psi)=\psi^{3}(1-\psi)^{3}$ and $\tilde{D}=\mathrm{Ra}/3$ for the case
of a slit Martin  _et al._ (2011). Equation (28) is solved numerically
following the method detailed by Martin et al. in Ref. Martin  _et al._
(2011). Figures 7(a–c) show this solution for the corresponding time
$t=1.6\times 10^{-2}$ superimposed with both the concentration field and the
velocity field. These data show a reasonable agreement, confirming the
negligible impact of diffusion on the solute transport.

To better describe this regime, Fig. 7(d) displays the average profiles
$\varphi_{0}(z,t)$ obtained from the 2D model at $\mathrm{Ra}=10^{5}$ and
$\mathrm{Sc}=10^{5}$ for time scales ranging from $t\simeq 0.003$ to $0.02$.
This plot shows that all the data almost collapse on a single curve when
plotted against $\xi=z/\sqrt{\tilde{D}t}$, which is correctly described by
$\psi(\xi)$ solution of Eq. (28). In this regime, one can again compute the
width of the mixing zone defined by Eq. (17) using $\psi(\xi)$, leading to:

$\displaystyle w\simeq 0.04879\sqrt{\mathrm{Ra}\,t}\,.$ (29)

This square-root spreading is plotted in Fig. 2(b) and accounts well for the
numerical data $w(t)$ obtained for $\mathrm{Ra}=10^{5}$ and time scales
ranging from $t\simeq 10^{-3}$ to $\simeq 10^{-1}$.

###  III.5 1D dispersion and late diffusion

Figure 8: (a) Concentration field $\varphi$, (b) component $u_{z}$ and (c)
component $u_{x}$ at $t=2.5$ ($\mathrm{Ra}=10^{5}$ and $\mathrm{Sc}=10^{5}$).
The black lines in (a) are isoconcentration lines. (d) Average concentration
profiles $\varphi_{0}(z,t)$ in the dispersion regime plotted against $\eta$
given by Eq. (34), $t$ ranges from $t=1$ to $2.5\times 10^{2}$ (10 curves).
The dashed line is Eq. (33). Inset: same profiles against $z$. (e) Rescaled
concentration profiles $\varphi_{0}(z,t)$ in the late diffusion regime,
$t=10^{4}$ up to $1.6\times 10^{5}$ (10 curves). The dashed line is Eq. (19).
Inset: same profiles against $z$.

For time scales $t\geq\mathcal{O}(1)$, 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=2.5$. In this regime, the transport is fully
controlled by the coupling between solute diffusion along the channel height
and buoyancy-driven advection along the channel main axis. This regime has
already been described in the literature since the pioneering work of Chatwin
and Erdogan, who studied the Taylor-Aris dispersion of a buoyant solute in a
pressure-driven flow Erdogan and Chatwin (1967), see also Smith (1976); Barton
(1976); Godfrey (1980); Maclean and Alboussière (2001); Salmon and Doumenc
(2020) and the review of Young and Jones on shear dispersion Young and Jones
(1991). In this regime, the extent $w$ of the concentration gradient along $z$
is large ($w\gg 1$), the buoyancy-driven flow is quasi-parallel [$u_{x}\ll
u_{z}$, see Figs. 8(b) and 8(c)], and the variations of the concentration
along $x$ are small [see the isoconcentration lines in Fig. 8(a)]. One can
therefore use the lubrication approximation to show that the density gradient
along $z$ drives a flow following:

$\displaystyle
u_{z}(x,z,t)=-\frac{\mathrm{Ra}}{12}\frac{\partial\varphi_{0}}{\partial
z}x(2x-1)(x-1)\,,$ (30)

see Appendix B. This flow adds a contribution to the dispersion which scales
as $\sim u_{z}^{2}$, as for the classical Taylor-Aris dispersion in a
Poiseuille flow. More rigorously, one can demonstrate that the average
concentration $\varphi_{0}(z,t)$ obeys the 1D dispersion equation:

$\displaystyle\frac{\partial\varphi_{0}}{\partial t}=\frac{\partial}{\partial
z}\left(D_{\mathrm{eff}}\frac{\partial\varphi_{0}}{\partial z}\right)\,,$ (31)

with:

$\displaystyle
D_{\mathrm{eff}}=1+\frac{1}{\alpha}\left(\text{Ra}\frac{\partial\varphi_{0}}{\partial
z}\right)^{2}\,,$ (32)

and $\alpha=362880$, see Appendix B. The non-linearity of the dispersive term
comes from the coupling between the concentration gradient and the flow,
unlike the case of the Taylor-Aris dispersion: strong gradient increases the
magnitude of the gravity current which in turn increases the dispersion of the
solute. We will return in detail to the derivation of these equations in the
next section when we tackle the 3D case of a microfluidic channel with a
rectangular cross-section.

Equation (31) can be made free of any parameter by defining
$t^{\star}=(\alpha/\mathrm{Ra}^{2})t$ and
$z^{\star}=(\sqrt{\alpha}/\mathrm{Ra})z$. It is then solved numerically with
the initial condition $\varphi_{0}(z,t=0)=\mathcal{H}(z)$ to compute the
solution for any Rayleigh number and $\alpha$ value. Figure 2(b) shows that
the width of the mixing zone $w(t)$ defined by Eq. (17) and computed from the
numerical resolution of Eq. (31), perfectly matches the data obtained from the
full 2D model for time scales $t\geq 0.1$. The component $u_{z}$ at $z=0$
computed from the 2D model is also very well-approximated by Eq. (30) using
the solution of the 1D dispersion model (data not shown).

Equation (31) along with the initial condition
$\varphi_{0}(z,t=0)=\mathcal{H}(z)$ has been studied by Maclean and
Alboussière Maclean and Alboussière (2001) who provided asymptotic
approximations of the solution. When buoyancy dominates the transport of the
solute, i.e. $D_{\mathrm{eff}}\gg 1$, Eq. (31) admits the self-similar
solution:

$\displaystyle\varphi_{0}=\frac{1}{2}+\frac{\eta}{2}\sqrt{\frac{1}{\sqrt{3}\pi}-\frac{\eta^{2}}{12}}+\frac{1}{\pi}\text{arcsin}\left(\frac{\sqrt{\pi}\eta}{2\times
3^{1/4}}\right)\,,$ (33)

with $\eta$ given by:

$\displaystyle\eta=\frac{\alpha^{1/4}}{\sqrt{{\rm Ra}}}\frac{z}{t^{1/4}}\,.$
(34)

Eq. (33) is valid for $\eta^{2}\leq 12/(\pi\sqrt{3})$ Maclean and Alboussière
(2001). Figure 8(d) displays this asymptotic self-similar solution for
$\mathrm{Ra}=10^{5}$ 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:

$\displaystyle
w\simeq\sqrt{\frac{\text{Ra}}{2\pi}}\left(\frac{3t}{\alpha}\right)^{1/4}\,,$
(35)

thus following $w\sim t^{1/4}$ as shown in Fig. 2(b).

At later times, the density gradient continuously decreases as the solute is
continuously dispersed along the channel, and diffusion dominates again the
transport of the solute, i.e. $D_{\mathrm{eff}}\simeq 1$. In this late
diffusion regime, we once again find a classical diffusion problem where
solutal free convection no longer plays a role, and the concentration profiles
are then given by Eq. (19). This is illustrated by Fig. 8(e) showing Eq. (19)
superimposed with the average concentration profiles computed from the 2D
model at long time scales. In this late regime, the width of the mixing zone
is again given by $w=\sqrt{t}$ as shown in Fig. 2(b), despite the buoyancy-
driven flow along the slit still given by Eq. (30). As explained in
Introduction, even if this flow has no effect on the solute gradient that
generates it, it still exists and may have an effect on less mobile species in
the case of complex liquid mixtures.

###  III.6 Dispersion diagram

I- early diffusion | $w\simeq\sqrt{t}$ | $\bar{u}_{z}(z=0,t)\simeq 0.00690\,\mathrm{Ra}$ | $t_{\mathrm{I}\to\mathrm{II}}\simeq 34550/\mathrm{Ra}^{2}$
---|---|---|---
II- early advection | $w\simeq 0.00538\,\text{Ra}\,t$ | $\bar{u}_{z}(z=0,t)\simeq 0.00690\,\mathrm{Ra}$ | $t_{\mathrm{II}\to\mathrm{III}}\simeq 82/\mathrm{Ra}$
III- late advection | $w\simeq 0.04879\sqrt{\mathrm{Ra}\,t}~{}$ | $\bar{u}_{z}(z=0,t)\simeq 0.03138\sqrt{\frac{\mathrm{Ra}}{t}}$ | $t_{\mathrm{III}\to\mathrm{IV}}\simeq 0.04$
IV- 1D dispersion | $w\simeq\sqrt{\frac{\text{Ra}}{2\pi}}\left(\frac{3t}{\alpha}\right)^{1/4}$ | $\bar{u}_{z}(z=0,t)\simeq\frac{\sqrt{\mathrm{Ra}}}{192\sqrt{\pi}}\left(\frac{\alpha}{3t}\right)^{1/4}$ | $t_{\mathrm{IV}\to\mathrm{V}}\simeq\frac{3}{(2\pi)^{2}\alpha}\,\mathrm{Ra}^{2}$
V- late diffusion | $w\simeq\sqrt{t}$ | $\bar{u}_{z}(z=0,t)\simeq\frac{1}{384\sqrt{\pi}}\,\ \frac{\mathrm{Ra}}{\sqrt{t}}$ |
Table 1: Transport regimes and transition times for a slit, $\alpha=362880$.

In the previous paragraphs, we identified a sequence of regimes of transport
of the solute, which are summarized in Table 1. For each regime, one can
compute the typical longitudinal velocity at $z=0$ defined by:

$\displaystyle\bar{u}_{z}(z=0,t)=-2\int_{0}^{1/2}u_{z}(x,z=0,t)\text{d}x,$
(36)

using in particular Eq. (22) for Regimes I and II, Ref. Martin  _et al._
(2011) for Regime III, and Eq. (30) along with Eq. (33) [resp. Eq. (19)] for
Regime IV (resp. Regime V). Results are displayed in the third column of Table
1.

The transition times between these regimes are estimated by matching the
different spreading laws $w$ vs. $t$ leading to the values provided in Table
1. These definitions lead to different scaling laws with $\mathrm{Ra}$ which
are also reported in Fig. 2. These transition times allow us to construct the
diagram presented in Fig. 9(a) in the plane $\mathrm{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 also coupled for
large $\mathrm{Ra}$ and small $\mathrm{Sc}$ to the momentum diffusion across
the slit, see Sec. III.2 and in particular Fig. 4(b).

Figure 9: (a) Diagram of the different regimes in the plane $\mathrm{Ra}$ vs.
$t$. The transition times $t$ vs. $\mathrm{Ra}$ are given in Table 1. The
vertical dotted lines correspond to the extent of the inertial regime $t\simeq
0.0571/\text{Sc}$ for $\mathrm{Sc}=10^{3}$ and $\mathrm{Sc}=10^{5}$, see Sec.
III.2. (b) Width of the mixing zone $w(t)$ for $\mathrm{Sc}=10^{5}$ and
$\mathrm{Ra}=10^{5}$ (red), $10^{4}$ (blue), $10^{3}$ (magenta). The black
dashed line is the purely diffusive spreading $w=\sqrt{t}$. The symbols are
the transition times given in Table 1 for $\mathrm{Ra}=10^{5}$ and
$\mathrm{Ra}=10^{4}$.

This diagram also shows that the effect of buoyancy vanishes for
$\mathrm{Ra}\lesssim 10^{3}$. To illustrate this point, Fig. 9(b) reports the
spreading laws $w$ vs. $t$ computed from the numerical resolution of the full
2D model for several Rayleigh numbers $\mathrm{Ra}=10^{5}$, $10^{4}$, $10^{3}$
and $\mathrm{Sc}=10^{5}$. These data clearly reveal that buoyancy has little
effect on the transport of the solute at all time scales for
$\mathrm{Ra}=10^{3}$. More quantitatively, the ratio $w/\sqrt{t}$ reaches a
maximum of only $\simeq 1.24$ at $t\simeq 0.1$ for $\mathrm{Ra}=10^{3}$. This
result answers the question initially asked in Introduction as it allows to
assess a numerical value to $\mathrm{Ra}$ corresponding to negligible buoyancy
in a microfluidic slit at all time scales, at least regarding the active
solute that generates the gravity current.

The regimes of early diffusion and early advection are only visible for high
Rayleigh numbers, $\mathrm{Ra}\geq 10^{4}$, and small time scales,
$t<10^{-3}$. For most microfluidic configurations investigating molecular
solutes, these regimes might be difficult to observe even in a thick slit, as
$t=10^{-3}$ in the above diagram does not exceed a few seconds for
$H=500~{}\mu$m and $D\geq 10^{-10}$ m2/s. On the other hand, both regimes of
late advection and 1D dispersion should be easily observed as the transition
time between these two regimes $t_{\mathrm{III}\to\mathrm{IV}}$ is a few
minutes for $D\simeq 10^{-10}$ m2/s and $H=500~{}\mu$m. In the case of
colloidal dispersions, the regimes of early diffusion and early advection
might be observable as much lower $D$ values lead to much longer time scales.
For instance, $t=10^{-3}$ is now a few minutes for $H=500~{}\mu$m and
$D=10^{-12}$ m2/s corresponding to colloids of radius $100$ nm dispersed in
water.

###  III.7 Analogy with the case of a 2D porous layer

Szulczewski and Juanes Szulczewski and Juanes (2013) studied in a different
context (geological sequestration of CO2 in an aquifer), a problem similar to
the one described in Fig. 1, but considering a vertically confined porous
layer of thickness $H$. Their theoretical model is also based on Eq. (7) to
describe the solute transport and Eq. (6) for the overall mass conservation,
but the pore velocity field $\mathbf{U}$ follows Darcy’s law given by:

$\displaystyle\mathbf{U}=-\frac{\kappa}{\rho_{0}\nu\epsilon}\left[\nabla
P-(\rho(\Phi)-\rho_{0})\mathbf{g}\right]\,,$ (37)

where $\kappa$ is the permeability of the permeable rock and $\epsilon$ its
porosity (see Sec. II.1 for the other notations). Notice that Navier-Stokes
equations Eq. (5) turns to Stokes equations when inertia is neglected, that
does not reduce to Darcy’s law Eq. (37). Indeed, Stokes equations include the
diffusion of the momentum over the scale $H$ of the microfluidic channel, that
brings a fundamental difference with Darcy’s law.

The characteritic velocity resulting from the Darcy’s law Eq. (37) reads

$\displaystyle U_{D}=\frac{\beta\Phi_{i}g\kappa}{\nu\epsilon}\,,$ (38)

and the corresponding Rayleigh number comparing diffusion and advection by the
gravity current is:

$\displaystyle\widetilde{\mathrm{Ra}}=\frac{\beta\Phi_{i}g\kappa
H}{\nu\epsilon D}\sim\frac{U_{D}H}{D}\,,$ (39)

thus highlighting a scaling law with the thickness different from the
microfluidic case, $\widetilde{\mathrm{Ra}}\propto H$ vs. $\mathrm{Ra}\propto
H^{3}$. Interestingly, the range of Rayleigh numbers $\widetilde{\mathrm{Ra}}$
involved in the context of CO2 sequestration Szulczewski and Juanes (2013)
still corresponds to the range of $\mathrm{Ra}$ studied in the present work.

Despite the differences between these two models, Szulczewski and Juanes also
reported five distinct transport regimes in the porous rock that have strong
similarities to those reported in Fig. 2(a) for a microfluidic slit, see their
names in Table 2. In particular, the concentration profiles depend on the same
self-similar variables in each regime, and the different transition times
(displayed in the last column of Table 1 for our model) obey the same scaling
laws with the Rayleigh number. We believe that these similarities are related
to the linearity of the Darcy and Stokes equations 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 the solute flux across the interface $z=0$ defined by:

$\displaystyle f(t)=-\int_{0}^{1}\left[\varphi
u_{z}-\left(\frac{\partial\varphi}{\partial
z}\right)\right]_{x,z=0}\text{d}x\,,$ (40)

as Szulczewski and Juanes also computed this quantity in each regime
Szulczewski and Juanes (2013). $f$ is estimated using Eqs. (19) in Regimes I
and V, Eq. (22) in Regime II, following Ref. Martin  _et al._ (2011) in Regime
III, and using Eq. (31) in Regime IV. As shown in Table 2, $f$ for the
microfluidic slit and $\tilde{f}$ for the porous layer show the same scaling
laws with $t$ and $\mathrm{Ra}$ (resp. $\widetilde{\mathrm{Ra}}$). With the
exception of the diffusion Regimes I and V, the large differences between the
numerical prefactors (up to two orders of magnitude) confirm that both
problems are fundamentally different. Relevant prefactors must be considered
in potential comparisons with experimental works.

Microfluidic slit, this work | 2D porous layer Szulczewski and Juanes (2013)
---|---
I- early diffusion | $f\simeq\frac{1}{2(\pi t)^{1/2}}$ | I- early diffusion | $\tilde{f}\simeq\frac{1}{2(\pi t)^{1/2}}$
II- early advection | $f\simeq 0.00345\,\mathrm{Ra}$ | II- S-slumping | $\tilde{f}\simeq 0.186\,\widetilde{\mathrm{Ra}}$
III- late advection | $f\simeq 0.0157\,\frac{\mathrm{Ra}^{1/2}}{t^{1/2}}$ | III- straight-line slumping | $\tilde{f}\simeq 0.125\,\frac{\widetilde{\text{Ra}}^{1/2}}{t^{1/2}}$
IV- 1D dispersion | $f\simeq 0.00321\,\frac{\mathrm{Ra}^{1/2}}{t^{3/4}}$ | IV- Taylor slumping | $\tilde{f}\simeq 0.0238\,\frac{\widetilde{\mathrm{Ra}}^{1/2}}{t^{3/4}}$
V- late diffusion | $f\simeq\frac{1}{2(\pi t)^{1/2}}$ | V- late diffusion | $\tilde{f}\simeq\frac{1}{2(\pi t)^{1/2}}$
Table 2: Solute flux defined by Eq. (40) for a microfluidic slit and in a 2D
porous layer Szulczewski and Juanes (2013).

##  IV The case of a 3D microfluidic channel

We now consider the case of a microfluidic channel with a square cross-
section, i.e. $\gamma=L/H=1$ in Fig. 1. Numerical simulations were performed
for three Rayleigh numbers $\mathrm{Ra}=10^{3}$, $10^{4}$, and $10^{5}$ with a
fixed Schmidt number $\mathrm{Sc}=10^{3}$.

###  IV.1 Transport regimes in the 3D case

Movie M2 supplied in the ESI shows the concentration and velocity fields in
the planes $y=0$ and $z=0$ obtained from the numerical resolution of the 3D
model for $\mathrm{Ra}=10^{5}$ and $\mathrm{Sc}=10^{3}$. This movie helps to
identify the same succession of regimes as for the slit case. From the full
numerical data, we computed again the width of the mixing zone $w(t)$ using
Eqs. (16) and (17) for three different Rayleigh numbers $\mathrm{Ra}=10^{3}$,
$10^{4}$ and $10^{5}$ and the same Schmidt number $\mathrm{Sc}=10^{3}$, see
Fig. 10. The dispersion curves obtained from the 3D model follow the same
trends as the 2D ones shown in Fig. 2(b). These data lead to the same
conclusion as for the slit case: buoyancy hardly affects the solute transport
at all time scales for $\mathrm{Ra}\leq 10^{3}$ (the maximal value of the
ratio $w/\sqrt{t}$ is only $\simeq 1.14$ at $t\simeq 0.1$ for
$\mathrm{Ra}=10^{3}$).

Figure 10: Width of the mixing zone $w(t)$ obtained in the case of a channel
with a square cross-section for $\mathrm{Sc}=10^{3}$ and $\mathrm{Ra}=10^{5}$
(red), $10^{4}$ (blue), $10^{3}$ (magenta). Asymptotic models for
$\mathrm{Ra}=10^{5}$: the yellow dashed line given by Eq. (41) corresponds to
the regime of late advection; the green dashed line is computed from the
numerical resolution of Eqs. (31) and (32) with $\alpha\simeq 739872$. The
black dashed line indicates the diffusion law $w=\sqrt{t}$.

These observations evidence that the description of the regimes presented
above applies again for the 3D case. Thereafter, we will not re-describe the
regimes of early diffusion and early advection, in particular because they are
short or even hardly observable in most microfluidic experimental
configurations. With respect to the late advection regime, the analysis
reported in Sec. III.4 can easily be adapted to the 3D case. The shape of the
pseudo-interface separating the solution and the solvent is computed using Eq.
(28) with $f(\psi)$ and $\tilde{D}$ corresponding to a channel with a square
cross-section, see Ref. Martin  _et al._ (2011) for details. We then
calculated the width $w(t)$ from the theoretical profiles, leading to:

$\displaystyle w\simeq 0.04104\sqrt{\mathrm{Ra}\,t}\,.$ (41)

This prediction fits well the data obtained from the 3D numerical simulation
at $\mathrm{Ra}=10^{5}$ for time scales ranging from $t\simeq 3\times 10^{-3}$
to $\simeq 3\times 10^{-2}$, see Fig. 10. The numerical prefactor in Eq. (41)
($0.04104$) is $16\%$ smaller than in Eq. (29) for the slit case ($0.04879$),
due to the increased viscous forces induced by the side walls in the 3D case.

At later time scales $t\geq\mathcal{O}(1)$, Movie M2 evidences that diffusion
almost homogenizes the concentration over the cross-section of the channel.
This is illustrated by Figs. 11(a) and 11(d) showing the concentration field
in the planes $y=0$ and $z=0$ at $t\simeq 4.7$. Figure 11, displaying also the
components of the velocity field in the same planes and at the same time,
evidences a quasi-parallel flow along $z$ (due to the large extent $w$ of the
longitudinal density gradient), and a secondary transverse flow as revealed by
the components $u_{x}$ and $u_{y}$ in the $z=0$ plane.

Figure 11: 3D case, $\mathrm{Ra}=10^{5}$, $\mathrm{Sc}=10^{3}$, and $t\simeq
4.7$. (a) Concentration field $\varphi(x,y=0,z,t)$. (b) Velocity field
$u_{z}(x,y=0,z,t)$ and (c) $u_{x}(x,y=0,z,t)$. (d) $\varphi(x,y,z=0,t)$, (e)
$u_{z}(x,y,z=0,t)$, (f) $u_{x}(x,y,z=0,t)$, and (g) $u_{y}(x,y,z=0,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 regime described 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 the literature for a
rectangular channel. It has been already calculated for a circular tube
Erdogan and Chatwin (1967) and a slit Young and Jones (1991), and we derive it
for a rectangular cross-section below in Sec. IV.2; (ii) As mentioned by
Chatwin and Ergogan Erdogan and Chatwin (1967) 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 in Sec. IV.3.

###  IV.2 1D dispersion regime in the 3D case

The fact that $w\gg 1$ calls for the ”Taylor-like” approach leading to the 1D
dispersion equation Eq. (31), see Appendix B for the slit case. In such a
theoretical approach Young and Jones (1991), transverse variations in
concentration are assumed to be small:

$\displaystyle\varphi(x,y,z,t)=\varphi_{0}(z,t)+\varphi_{1}(x,y,z,t)\,,$ (42)

with $\varphi_{0}$ given by Eq. (16) and $\varphi_{1}\ll\varphi_{0}$.
Averaging the transport equation Eq. (12) over the cross-section of the
channel leads to:

$\displaystyle\frac{\partial\varphi_{0}}{\partial
t}+\frac{\partial<u_{z}\varphi_{1}>}{\partial
z}=\frac{\partial^{2}\varphi_{0}}{\partial z^{2}}\,,$ (43)

similarly to Eq. (65) for the case of a slit. Subtracting this relation to Eq.
(12) gives:

$\displaystyle\frac{\partial\varphi_{1}}{\partial
t}+u_{z}\frac{\partial\varphi_{0}}{\partial
z}+\mathbf{u}.\nabla\varphi_{1}-\frac{\partial<u_{z}\varphi_{1}>}{\partial
z}=\frac{\partial^{2}\varphi_{1}}{\partial
x^{2}}+\frac{\partial^{2}\varphi_{1}}{\partial
y^{2}}+\frac{\partial^{2}\varphi_{1}}{\partial z^{2}}\,,$ (44)

as Eq. (66). For $t\gg 1$ and $w\gg 1$, Eq. (44) yields at leading order:

$\displaystyle u_{z}\frac{\partial\varphi_{0}}{\partial
z}+u_{y}\frac{\partial\varphi_{1}}{\partial
y}+u_{x}\frac{\partial\varphi_{1}}{\partial
x}\simeq\frac{\partial^{2}\varphi_{1}}{\partial
x^{2}}+\frac{\partial^{2}\varphi_{1}}{\partial y^{2}}\,.$ (45)

In the 2D case of a slit, $u_{y}=0$ and the continuity equation Eq. (11)
imposes $u_{x}\sim u_{z}/w\ll u_{z}$. The term $u_{x}\partial_{x}\varphi_{1}$
in the above relation is thus negligible, and the derivation of the dispersion
equation Eq. (31) is straightforward, see Appendix B. In the 3D case
nevertheless, the continuity equation does 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 $\partial_{x}u_{x}+\partial_{y}u_{y}=0$.
It is therefore a priori not possible to derive the dispersion equation Eq.
(31) unless one neglects the terms $u_{x}\partial_{x}\varphi_{1}$ and
$u_{y}\partial_{y}\varphi_{1}$ in Eq. (44). This point had been mentioned by
Chatwin and Erdogan who studied the classical Taylor-Aris problem of the
dispersion of a buoyant solute flowing in a circular tube Erdogan and Chatwin
(1967). They even showed that the lateral mixing of the solute due to these
secondary transverse flows could lead to a decrease of the overall solute
dispersion in a pressure-driven flow for some range of the Rayleigh number,
see also Refs. Barton (1976); Smith (1976).

The data shown in Fig. 11 along with Movie M2 in the ESI evidence that the
maximal magnitude of these secondary flows is about $u_{x}\simeq u_{y}\simeq
1$ for $\mathrm{Ra}=10^{5}$, and we will assume as a first step that they do
not significantly affect the transport of the solute. This assumption allows
us to neglect the terms $u_{y}\partial_{y}\varphi_{1}$ and
$u_{x}\partial_{x}\varphi_{1}$ in Eq. (44) as compared to the diffusive terms
$\partial_{x}^{2}\varphi_{1}$ and $\partial_{y}^{2}\varphi_{1}$, leading to:

$\displaystyle u_{z}\frac{\partial\varphi_{0}}{\partial
z}\simeq\frac{\partial^{2}\varphi_{1}}{\partial
x^{2}}+\frac{\partial^{2}\varphi_{1}}{\partial y^{2}}\,,$ (46)

similarly to Eq. (67) for the slit.

To compute the longitudinal flow $u_{z}$, we first assume that it is described
by the lubrication approximation because the extent of the mixing zone is
large ($w\gg 1$):

$\displaystyle\frac{\partial^{2}u_{z}}{\partial
x^{2}}+\frac{\partial^{2}u_{z}}{\partial y^{2}}=\frac{\partial
p_{\ell}}{\partial z}\,,$ (47) $\displaystyle 0=\frac{\partial
p_{\ell}}{\partial y}\,,$ (48) $\displaystyle 0=\frac{\partial
p_{\ell}}{\partial x}+\mathrm{Ra}\varphi_{0}\,,$ (49)

where $p_{\ell}$ 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:

$\displaystyle\frac{\partial^{2}u_{z}}{\partial
x^{2}}+\frac{\partial^{2}u_{z}}{\partial
y^{2}}=\mathrm{Ra}\frac{\partial\varphi_{0}}{\partial
z}\left(\frac{1}{2}-x\right)\,,$ (50)

the term $1/2$ ensuring that $<u_{z}>=0$.

From Eqs. (46) and (50), it is therefore possible to calculate both $u_{z}$
and $\varphi_{1}$, and finally to compute the dispersion term
$<u_{z}\varphi_{1}>$ in Eq. (43) leading ultimately to the same dispersion
equation Eq. (31) but with a different $\alpha$ value for $D_{\mathrm{eff}}$
in Eq. (32). Appendix C reports solutions of Eqs. (46) and (50) using Fourier
series for a rectangular channel of arbitrary aspect ratio $\gamma$, see Eqs.
(70) and (71). These calculations allow us to compute the numerical values of
$\alpha$ as a function of $\gamma$, see Fig. 15 in Appendix C. For a square
cross-section, $\alpha\simeq 739872$, and $\alpha$ tends towards the value
derived for the slit $\alpha=362880$ for infinitely thin channel, unlike the
Taylor-Aris coefficient in a rectangular channel Chatwin and Sullivan (1982);
Ajdari  _et al._ (2006).

Figures 12(a) and 12(b) shows both the component $u_{z}$ and $\varphi_{1}$ in
the plane $z=0$ calculated using Eqs. (70) and (71) and the 1D dispersion
model Eq. (31) at $t\simeq 4.7$. These data well match the numerical data
reported in Fig. 11(d) and (e) at the same time, see also the comparisons
along the line $y=z=0$ shown in Figs. 12(c) and (d). Consequently, the
section-averaged concentration profiles $\varphi_{0}(z,t)$ are correctly
predicted by the 1D model (data not shown). These profiles make it possible to
compute the width of the mixing zone $w(t)$ [Eq. (17)] plotted in Fig. 10
along with the data obtained from the 3D numerical simulation for
$\mathrm{Ra}=10^{5}$. The 1D dispersion approach accounts well for the
spreading of the solute for $t\geq 0.1$. Moreover, Eq. (35) giving the $w\sim
t^{1/4}$ behavior, allows us to compare the 2D and 3D configurations in this
regime. Eq. (35) indeed indicates that the ratio
$w_{\text{2D}}(t)/w_{\text{3D}}(t)$ between the 2D and the 3D cases is
$(\alpha_{\text{2D}}/\alpha_{\text{3D}})^{1/4}\simeq(362880/739872)^{1/4}\simeq
0.84\leq 1$, again due to the side walls increasing viscous forces in the 3D
case.

Figure 12: (a) $u_{z}(x,y,z=0)$ and (b) $\varphi_{1}(x,y,z=0)$ calculated
using Eqs. (70) and (71) at $t\simeq 4.7$, see the corresponding numerical
data shown in Fig. 11(d) and (e). The longitudinal density gradient is
calculated by the numerical resolution of the dispersion equation Eq. (31).
(c) and (d) show the comparisons along the line $y=z=0$ between the numerical
data reported in Fig. 11(d) and (e) (squares) and the theoretical profiles
given by Eqs. (70) and (71) (black lines).

###  IV.3 Transverse flows and validity range of the 1D dispersion equation

The comparisons shown in Sec. IV.2 evidence the validity of the 1D dispersion
equation to describe the overall solute transport at least for
$\mathrm{Ra}\lesssim 10^{5}$, although the transverse flow is not considered
in the 1D model. In the following, we derive the expression of the transverse
flow with the assumption that it does not significantly affect the
concentration field. Then, we determine the critical Rayleigh number for which
this hypothesis no longer holds.

Our calculations make it possible to compute the transverse variations in
concentration $\varphi_{1}$, see Fig. 12(b). These variations evidence density
gradients along $y$ due to the presence of the lateral walls which impact the
longitudinal flow, and thus the solute distribution. These density gradients
are responsible for a secondary flow because they cause a pressure gradient
along $y$. We compute this flow assuming that it is locally invariant along
$z$ because of the large extent of the mixing zone ($w\gg 1$) and solenoidal,
i.e. $\partial_{y}u_{y}+\partial_{x}u_{x}=0$, to ensure the global mass
conservation. This flow is thus solution of:

$\displaystyle 0=\frac{\partial p_{t}}{\partial z}\,,$ (51)
$\displaystyle\frac{\partial^{2}u_{y}}{\partial
x^{2}}+\frac{\partial^{2}u_{y}}{\partial y^{2}}=\frac{\partial p_{t}}{\partial
y}\,,$ (52) $\displaystyle\frac{\partial^{2}u_{x}}{\partial
x^{2}}+\frac{\partial^{2}u_{x}}{\partial y^{2}}=\frac{\partial p_{t}}{\partial
x}+\mathrm{Ra}\varphi_{1}\,,$ (53)

where $p_{t}$ is the pressure field associated to this 2D solenoidal flow. As
in the work of Chatwin and Erdogan Erdogan and Chatwin (1967), the global flow
is therefore assumed to be the superposition of the longitudinal flow
described by Eqs. (47-49) with the secondary transverse flow determined by
Eqs. (51-53).

Figure 13: Components $u_{x}$ (a) and $u_{y}$ (b) scaled by
$10^{-8}\,\left(\mathrm{Ra}\frac{\partial\varphi_{0}}{\partial z}\right)^{2}$
computed from the resolution of the biharmonic equation Eq (76). (c)
Corresponding transverse flow field in the plane $x$-$y$. This solenoidal flow
field results from the density gradients given by $\varphi_{1}$, see Eq. (76).

An analytical approximation of the solution of Eqs. (51-53) is computed in
Appendix D for a rectangular cross-section of arbitrary aspect ratio $\gamma$
using a method described by Shankar et al. Shankar (1993); Shankar  _et al._
(2002). Figure 13 shows the theoretical prediction of this secondary
transverse flow for $\gamma=1$. The vector velocity field shown in Fig. 13(c)
simply corresponds to the solenoidal flow associated to the density gradients
revealed in Fig. 12(b). The predicted flow pattern compares well with the
numerical data reported in Fig. 11(f) and (g), and our calculations show in
particular that the maximal values of the components $u_{x}$ and $u_{y}$ are:

$\displaystyle\mathrm{max}(u_{x})\simeq\mathrm{max}(u_{y})\simeq 14\times
10^{-8}\left(\mathrm{Ra}\frac{\partial\varphi_{0}}{\partial z}\right)^{2}\,.$
(54)

To go a step further into the comparison, Fig. 14 displays the maximal values
of the components $u_{x}$, $u_{y}$ and $u_{z}$ in the transverse plane $z=0$
vs. time $t$ obtained from the 3D model for three Rayleigh numbers,
$\mathrm{Ra}=10^{3}$, $10^{4}$, and $10^{5}$.

Figure 14: Maximal value of the components $u_{z}$ (red), $u_{x}$ (blue) and
$u_{y}$ (magenta) in the plane $z=0$ computed from the numerical 3D model at
(a) $\mathrm{Ra}=10^{3}$, (b) $10^{4}$, and (c) $10^{5}$ ($\mathrm{Sc}=10^{3}$
for these three cases). The black lines are the prediction given by Eq. (70),
and the dotted line the prediction given by Eq. (54). In these theoretical
predictions, the longitudinal density gradient $\partial_{z}\varphi_{0}$ at
$z=0$ is calculated from the numerical solution of the 1D dispersion Eq. (31).

The theoretical prediction given by Eq. (70) for the longitudinal flow $u_{z}$
correctly fits the numerical data even at $\mathrm{Ra}=10^{5}$ for $t\geq 1$.
The theoretical transverse components given by Eq. (54) accounts well for the
data at $\mathrm{Ra}=10^{3}$ and $\mathrm{Ra}=10^{4}$ [with the gradient of
$\varphi_{0}$ at $z=0$ computed from the numerical solution of Eq. (31)], but
significant discrepancies are observed for $\mathrm{Ra}=10^{5}$. For this
Rayleigh number, $u_{x}\simeq u_{y}\simeq 1$ when entering the 1D dispersion
regime ($t\geq 1$), and advection by the secondary transverse flow is no more
negligible compared to the transverse diffusion. It is remarkable to see that
the theoretical prediction given by Eq. (54) makes it possible to predict this
transition quantitatively. More precisely, the concentration profiles in the
dispersion regime are correctly described by the asymptotic approximation Eq.
(33) when buoyancy dominates, and one can thus compute the longitudinal
gradient at $z=0$ to write Eq. (54) as follows:

$\displaystyle\mathrm{max}(u_{x})\simeq\mathrm{max}(u_{y})\simeq 2.2\times
10^{-5}\,\frac{\mathrm{Ra}}{\sqrt{t}}\,,$ (55)

for a square cross-section ($\alpha\simeq 739872$). Therefore, the impact of
transverse flows on the solute transport is expected to be negligible for
$\mathrm{max}(u_{x})\simeq\mathrm{max}(u_{y})\leq 1$ leading to the criterion:

$\displaystyle t\geq 4.9\times 10^{-10}\,\mathrm{Ra}^{2}\,,$ (56)

about $t\geq 4.9$ for $\mathrm{Ra}=10^{5}$, in a remarkable agreement with the
data of Fig. 14(c) showing a discrepancy between the numerical solutions
($u_{x}$,$u_{y}$) and the predictions given by Eq. (54) for $t\leq 10$.
Nevertheless, the transverse flow does not yet significantly change the
longitudinal dispersion even for $\mathrm{Ra}=10^{5}$ at $t>1$ (see the good
agreement between the 1D dispersion model and the 3D numerical simulation
noted in Sec. IV.2). However, our theoretical predictions suggest that a
significant impact on the overall solute dispersion is expected for higher
Rayleigh numbers. The secondary transverse flow probably leads in this case to
a decrease of the overall solute dispersion because it contributes to mix the
solute laterally, as also shown by Chatwin and Erdogan in a different context,
the Taylor-Aris dispersion of a buoyant solute in a pressure-driven flow
Erdogan and Chatwin (1967); Barton (1976); Smith (1976). In most microfluidic
experimental configurations, these effects can a priori be neglected because
$\mathrm{Ra}\leq 10^{5}$.

##  V Conclusions

In the present work, we have studied in detail the impact of buoyancy on
solute spreading in two 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 the transport
equations. One of the main results of our study is to show that for
$\mathrm{Ra}\leq 10^{3}$, solutal free convection does not impact solute
diffusion at all time scales. Beyond this result, our theoretical predictions
give also for larger Rayleigh numbers, the time scales (or density gradients)
for which buoyancy no longer impacts molecular diffusion, see the diagram in
Fig. 9. Moreover, these same theoretical predictions allow to estimate
analytically the gravity currents, whatever their role on solute transport. It
is worth remembering that these flows may impact the transport of other
species dispersed 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=100~{}\mu$m, free convection is not expected to impact the diffusive mixing
because $\mathrm{Ra}\simeq 230$. Nevertheless, the typical longitudinal
velocity defined by Eq. (36) is about $\bar{u}_{z}\simeq 25~{}\mu$m/s in the
early regimes of diffusion and advection, and still $\bar{u}_{z}\simeq
5~{}\mu$m/s for the time scale $T=H^{2}/D$ (see Table I), and could
significantly advect less mobile species dispersed in the solutions. The
precise control of transport conditions in microfluidic geometries thus
possibly opens the way to flow control induced by solute gradients.

For the 3D case of a rectangular cross-section channel, our work brings for
the first time (to our knowledge) estimates of the 1D dispersion coefficient
describing the transport of the solute at long time scales for any aspect
ratio $\gamma=L/H$. Our work also highlights a subtle point related to 3D
geometries: the order of magnitude of the transverse flows cannot be
determined from the lubrication approximation alone. In the case studied here,
these flows, induced by transverse density gradients, remain moderate up to
$\mathrm{Ra}=10^{5}$, and the overall solute spreading is correctly described
by a 1D dispersion equation. However, our work predicts that these flows could
play a role at higher Rayleigh numbers for experimental situations outside the
field of application of microfluidics.

It could also be relevant to study more in details the case of shallow
channels commonly encountered in microfluidic applications, and in particular
to study more finely the transition between a channel with a large aspect
ratio, $\gamma\gg 1$, and the slit. Indeed, the case $\gamma\gg 1$ deserves
more attention because a new time scale appears, the diffusion time over the
width of the channel $\sim\gamma^{2}$, see Ref. Ajdari  _et al._ (2006)
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 the Boussinesq approximation, see Eqs. (5–7). Microfluidic
technologies allow a very fine control of the transport conditions (especially
mass and momentum), and thus to study interdiffusion in more complex mixtures.
It would then be useful to go beyond the model described by Eqs. (5–7) to
include this complexity: change in viscosity and diffusion coefficient as a
function of concentration, role of the reference frame (volume velocity / mass
velocity) Joseph  _et al._ (1996); Brenner (2005), other transport mechanisms
(e.g. diffusio-osmosis), etc. In this context, we hope that our work will make
it possible to disentangle the role played by solutal free convection from
other transport phenomena.

## Acknowledgements

This work was performed using HPC resources from the ”Mésocentre” computing
center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported
by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/). JBS
also thanks Y. Hallez for discussions concerning the lock-exchange problem and
ANR OSMOCHIP (ANR-18-CE06-0021) as well as Solvay 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 velocity field is solution of the Stokes equation Eq. (20).
Introducing the stream function $\psi(x,z)$ defined by:

$\displaystyle u_{x}=\frac{\partial\psi}{\partial
z}~{}~{}\mathrm{and}~{}~{}u_{z}=-\frac{\partial\psi}{\partial x}\,,$ (57)

Eq. (20) is equivalent to:

$\displaystyle\Delta^{2}\psi=\mathrm{Ra}\delta(z)\,,$ (58)

with $\delta(z)$ the Dirac function and $\Delta^{2}$ the biharmonic operator
in 2D. The no-slip boundary conditions impose:

$\displaystyle\frac{\partial\psi}{\partial
z}(x=0~{}\mathrm{and}~{}1,z)=0~{}~{}\mathrm{and}~{}~{}\frac{\partial\psi}{\partial
x}(x=0~{}\mathrm{and}~{}1,z)=0\,.$ (59)

The flow is expected to vanish far from $z=0$ and we impose
$\psi(x,z\to\pm\infty)=0$. Because
$\text{d}\psi=-u_{z}\text{d}x+u_{x}\text{d}z$, integrating $\text{d}\psi$ over
a contour along $x=0$ or $x=1$ starting from $z\to\pm\infty$ imposes:

$\displaystyle\psi(x=0~{}\mathrm{and}~{}1,z)=0\,.$ (60)

We define the Fourier transform of $\psi(x,z)$ by:

$\displaystyle\tilde{\psi}(x,k)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}{\rm
d}z\psi(x,z)e^{ikz}\,.$ (61)

Eq. (58) turns to:

$\displaystyle\frac{\partial^{4}\tilde{\psi}}{\partial
x^{4}}-2k^{2}\frac{\partial^{2}\tilde{\psi}}{\partial
x^{2}}+k^{4}\tilde{\psi}=\frac{\mathrm{Ra}}{\sqrt{2\pi}}\,.$ (62)

The solution of the ordinary differential equation is:

$\displaystyle\tilde{\psi}(x,k)=\mathrm{Ra}\frac{k+k(x-1)\cosh(kx)-kx\cosh(k-kx)+\sinh(k)-\sinh(kx)-\sinh(k-kx)}{\sqrt{2\pi}\,k^{4}(k+\sinh(k))},$
(63)

for the above boundary conditions. The stream function is then computed from
the inverse Fourier transform of $\tilde{\psi}(x,k)$ and the velocity field is
found using Eqs. (57). These calculations lead in particular to Eqs. (21) and
(22) given in Sec. III.3. Due to the symmetry of the equations along the
$x=1/2$ plane, the maximum of $u_{x}$ occurs at $x=1/2$ for all $z$ values.

##  Appendix B Derivation of the advection-dispersion equation for the case of
a slit

We assume:

$\displaystyle\varphi(x,z,t)=\varphi_{0}(z,t)+\varphi_{1}(x,z,t)\,,$ (64)

with $\varphi_{0}$ given by Eq. (16) and $\varphi_{1}\ll\varphi_{0}$.
Averaging the transport equation Eq. (12) over the height of the slit with the
help of the continuity relation Eq. (11), leads to:

$\displaystyle\frac{\partial\varphi_{0}}{\partial
t}+\frac{\partial<u_{z}\varphi_{1}>}{\partial
z}=\frac{\partial^{2}\varphi_{0}}{\partial z^{2}}\,.$ (65)

Subtracting this relation to Eq. (12) yields:

$\displaystyle\frac{\partial\varphi_{1}}{\partial
t}+u_{z}\frac{\partial\varphi_{0}}{\partial
z}+\mathbf{u}.\nabla\varphi_{1}-\frac{\partial<u_{z}\varphi_{1}>}{\partial
z}=\frac{\partial^{2}\varphi_{1}}{\partial
x^{2}}+\frac{\partial^{2}\varphi_{1}}{\partial z^{2}}\,.$ (66)

If we now assume that the scale $w$ of the gradient along $z$ verifies $w\gg
1$, the continuity equation imposes $u_{x}\sim u_{z}/w$, and Eq. (66) reduces
to:

$\displaystyle u_{z}\frac{\partial\varphi_{0}}{\partial
z}\simeq\frac{\partial^{2}\varphi_{1}}{\partial x^{2}}\,,$ (67)

provided that $t\gg 1$ and $\varphi_{1}\ll\varphi_{0}$. Similarly, by
neglecting the inertial term in Eq. (10) and assuming again $w\gg 1$, one find
Eq. (30) for the component $u_{z}$, and $u_{x}$ given by the continuity
equation Eq. (11). Integration of Eq. (67) along with the impermeability
boundary condition and $<\varphi_{1}>=0$ gives:

$\displaystyle\varphi_{1}=-\frac{\text{Ra}}{1440}\left(\frac{\partial\varphi_{0}}{\partial
z}\right)^{2}\left(12x^{5}-30x^{4}+20x^{3}-1\right)\,.$ (68)

Inserting Eq. (68) into Eq. (65) leads to the dispersion equation Eq. (31)
with the dispersion coefficient given by Eq. (32). Comparison of the solution
of Eq. (31) with the data obtained from the full 2D model shows that the 1D
dispersion model is actually valid as soon as $t\geq 1$.

##  Appendix C Derivation of the advection-dispersion equation for the case of
a microfluidic channel with a rectangular cross-section

The solution of Eq. (50) is found by a Fourier sum noticing first that:

$\displaystyle\frac{1}{2}-x=\sum_{n=2,\mathrm{even}}^{\infty}\frac{2}{n\pi}\sin(n\pi
x)\,~{}~{}\text{for}~{}~{}x\in]0-1[\,.$ (69)

The Fourier series representing the solution of Eq. (50) along with the no-
slip boundary conditions at the solid walls is then:

$\displaystyle u_{z}=\mathrm{Ra}\frac{\partial\varphi_{0}}{\partial
z}\sum_{n=2,\mathrm{even}}^{\infty}\frac{2}{(n\pi)^{3}}\left(\frac{\cosh(\pi
ny)}{\cosh(\pi n\gamma/2)}-1\right)\sin(n\pi x)\,,$ (70)

with $\gamma=L/H$, see Ref. Bruus (2007) for a similar calculation of the
pressure-driven velocity profile in a rectangular channel. The Fourier series
representing $\varphi_{1}$ can now be found using Eq. (46) along with the
impermeability boundary conditions at the solid walls and the constraint
$<\varphi_{1}>=0$, leading after calculations to:

$\displaystyle\varphi_{1}=\frac{8\mathrm{Ra}}{\pi^{6}}\left(\frac{\partial\varphi_{0}}{\partial
z}\right)^{2}\sum_{n=1,\mathrm{odd}}^{\infty}\sum_{p=2,\mathrm{even}}^{\infty}\frac{\cos(n\pi
x)}{p^{2}(p^{2}-n^{2})^{2}}$ $\displaystyle\bigg{[}\frac{\cosh(\pi
py)}{\cosh(\pi p\gamma/2)}$ (71) $\displaystyle-\frac{p}{n}\tanh(\pi
p\frac{\gamma}{2})\frac{\cosh(\pi ny)}{\sinh(\pi
n\gamma/2)}+\frac{p^{2}-n^{2}}{n^{2}}\bigg{]}.$

The dispersion term $<u_{z}\varphi_{1}>$ in Eq. (43) can now be evaluated
leading finally to the dispersion equation Eq. (31) with $D_{\mathrm{eff}}$
given by Eq. (32). The numerical prefactor $\alpha$ is given by:

$\displaystyle\frac{1}{\alpha}=\frac{1}{362880}-\frac{64}{\gamma\pi^{11}}\sum_{n=1,\mathrm{odd}}^{\infty}\sum_{p=2,\mathrm{even}}^{\infty}\sum_{q=2,\mathrm{even}}^{\infty}\frac{1}{n^{2}p^{2}q^{2}(p^{2}-n^{2})^{2}(q^{2}-n^{2})}\bigg{[}\frac{n^{2}\left(p\tanh\left(\frac{\pi\gamma
p}{2}\right)-q\tanh\left(\frac{\pi\gamma q}{2}\right)\right)}{p^{2}-q^{2}}$
$\displaystyle-\frac{np\tanh\left(\frac{\pi\gamma
p}{2}\right)\left(n-q\coth\left(\frac{\pi\gamma
n}{2}\right)\tanh\left(\frac{\pi\gamma
q}{2}\right)\right)}{n^{2}-q^{2}}+\frac{(p^{2}-n^{2})\tanh\left(\frac{\pi\gamma
q}{2}\right)}{q}+\frac{(p^{2}-n^{2})\tanh\left(\frac{\pi\gamma
p}{2}\right)}{p}\bigg{]}\,,$ (72)

and therefore only depends on the aspect ratio $\gamma$ of the channel.

Asymptotic approximations of $\alpha$ can be found for $\gamma\to 0$:

$\displaystyle\frac{1}{\alpha}=\frac{\gamma^{4}}{17280}+\mathcal{O}(\gamma^{6})\,.$
(73)

For a wide slit, i.e. $\gamma\gg 1$, the terms $\tanh$ and $\coth$ in Eq. (72)
are close to $\simeq 1$, and one thus find:

$\displaystyle\frac{1}{\alpha}\simeq\frac{1}{362880}-\frac{1}{681432\gamma}\,,$
(74)

where the factor $\simeq 681432$ is found by calculating numerically the sum
in Eq. (72). Figure 15 displays the values of $\alpha$ calculated using Eq.
(72) for several aspect ratios $\gamma$ along with the asymptotic behaviors
Eqs. (73) and (74).

Figure 15: Prefactor $\alpha$ in Eq. (32) for a rectangular cross-section of
aspect ratio $\gamma$. The dashed dotted line is the asymptotic bahavior Eq.
(73) valid for $\gamma\to 0$. The magenta line is the approximation Eq. (74)
valid for $\gamma\gg 1$. The dashed line is the value for a slit
$\alpha=362880$. Inset: zoom on the $\gamma$ range $1$–$4$.

It should be noted that the approximation Eq. (74) valid for thin channels
$\gamma\gg 1$ yields correct estimates of $\alpha$ within $<0.5\%$ for
$\gamma\geq 1.5$. The value for a square cross-section is $\alpha\simeq
739872$.

##  Appendix D Derivation of the secondary transverse flow

We introduce the stream function $\psi(x,y)$ defined by:

$\displaystyle u_{x}=\frac{\partial\psi}{\partial
y}~{}~{}\mathrm{and}~{}~{}u_{y}=-\frac{\partial\psi}{\partial x}\,.$ (75)

Eqs. (52) and (53) turn to the inhomogeneous biharmonic equation:

$\displaystyle\left(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial}{\partial
y^{2}}\right)^{2}\psi=\mathrm{Ra}\frac{\partial\varphi_{1}}{\partial y}\,,$
(76)

with $\varphi_{1}$ given by Eq. (71). The no-slip boundary conditions at the
solid walls impose:

$\displaystyle\frac{\partial\psi}{\partial
y}(x,y=\pm\gamma/2)=\frac{\partial\psi}{\partial
x}(x=0~{}\mathrm{and}~{}1,y)=\psi(x,y=\pm\gamma/2)=\psi(x=0~{}\mathrm{and}~{}1,y)=0\,.$
(77)

It can be seen that the same equations govern viscous flows induced by
inhomogeneous temperature fields in a rectangular container, and we will use
the method described in Refs. Shankar (1993); Shankar  _et al._ (2002) to
estimate the solution of Eq. (76). In brief, the general solution of Eq. (76)
is written as $\psi=\psi_{\mathrm{in}}+\psi_{\mathrm{in}}$ with
$\psi_{\mathrm{in}}$ any solution of the inhomogeneous equation (i.e. with the
right-hand term), and $\psi_{\mathrm{ho}}$ the solution of the homogenous
biharmonic equation, so that $\psi$ fulfills the above boundary conditions.
$\psi_{\mathrm{in}}$ can be found simply by a Fourier series expansion
following $\varphi_{1}$, and $\psi_{\mathrm{ho}}$ can be found using a direct
eigenfunction expansion, see below and Refs. Shankar (1993); Shankar  _et al._
(2002) for details.

For simplicity, we re-write $\varphi_{1}$ as:

$\displaystyle\varphi_{1}=\mathrm{Ra}\left(\frac{\partial\varphi_{0}}{\partial
z}\right)^{2}\sum_{n=1,\mathrm{odd}}^{\infty}\sum_{p=2,\mathrm{even}}^{\infty}h_{n,p}(y)\cos(n\pi
x)\,,$ (78)

with $h_{n,p}(y)$ given in Eq. (71). An inhomogeneous solution
$\psi_{\mathrm{in}}$ of Eq. (76) can be easily found using the following
Fourier series:

$\displaystyle\psi_{\mathrm{in}}=\left(\mathrm{Ra}\frac{\partial\varphi_{0}}{\partial
z}\right)^{2}\sum_{n=1,\mathrm{odd}}^{\infty}\sum_{p=2,\mathrm{even}}^{\infty}a_{n,p}(y)\cos(n\pi
x)\,,$ (79)

with $a_{n,p}(y)$ solutions of the following ordinary differential equations:

$\displaystyle
a^{(4)}_{n,p}(y)-2\pi^{2}n^{2}a^{(2)}_{n,p}(y)+\pi^{4}n^{4}a_{n,p}(y)=h^{(1)}_{n,p}(y)\,.$
(80)

The solutions of these equations with the following boundary conditions:

$\displaystyle a^{(1)}_{n,p}(\pm\gamma/2)=a_{n,p}(\pm\gamma/2)=0\,,$ (81)

lead to a Fourier series representing $\psi_{\mathrm{in}}$. For the sake of
clarity, the functions $a_{n,p}(y)$ are not written here. This Fourier series
fullfills all the boundary conditions given by Eq. (77) except
$\psi_{\mathrm{in}}(x=0,y)=-\psi_{\mathrm{in}}(x=1,y)=\kappa(y)\neq 0$.

The next step consists therefore in finding the solution $\psi_{\mathrm{ho}}$
of the homogeneous biharmonic equation with the boundary conditions given by
Eq. (77) but with
$\psi_{\mathrm{ho}}(x=0,y)=-\psi_{\mathrm{ho}}(x=1,y)=-\kappa(y)$. As the
symmetry of the problem imposes that $\psi$ is an odd function of $y$, we will
use the following odd eigunfunctions $\exp(\pm\lambda_{n}x)\phi_{n}(y)$ with:

$\displaystyle\phi_{n}(y)=y\cos(\lambda_{n}y)-\frac{\gamma}{2}\cot\left(\frac{\gamma\lambda_{n}}{2}\right)\sin(\lambda_{n}y)\,,$
(82)

and $\lambda_{n}$ the complex roots of the transcendental equation:

$\displaystyle\sin(\lambda_{n})=\lambda_{n}\,,$ (83)

which can be estimated by the Newton’s method. These eigenfunctions verify the
boundary conditions expected at $y=\pm\gamma/2$. As $\psi$ is a real and odd
function of $x$ with respect to $x=1/2$, an eigenfunction expansion for
$\psi_{\mathrm{ho}}(x,y)$ is:

$\displaystyle\psi_{\mathrm{ho}}(x,y)=\left(\mathrm{Ra}\frac{\partial\varphi_{0}}{\partial
z}\right)^{2}\sum_{n=1}^{\infty}\mathcal{R}e\left[b_{n}\phi_{n}(y)\left(e^{-\lambda_{n}x}-e^{-\lambda_{n}(1-x)}\right)\right]$
(84)

where the complex numbers $b_{n}$ have to be determined from the boundary
conditions at $x=0$ and $x=1$. To get an approximate solution, we proceed as
proposed by Shankar in Ref. Shankar (1993). First, the sum in Eq. (84) is
truncated to the first $N$ terms. Then, a least-squares procedure is used to
find the coefficients $b_{n}$ which yield the best expected boundary
conditions at $x=0$ for $m$ equidistant points over the interval
$y=[0-\gamma/2]$. This procedure corresponding to the resolution of $2N$
linear algebraic equations is performed using Mathematica. The coefficients
$b_{n}$ are rapidly converging and only a few eigenfunctions are needed to get
an accurate estimate of $\psi_{\mathrm{ho}}(x,y)$.

Figures 13(a) and 13(b) display the components $u_{y}$ and $u_{x}$ computed
from the stream function $\psi=\psi_{\mathrm{in}}+\psi_{\mathrm{ho}}$
calculated using the above procedure for $\gamma=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 the gravity along $x$, and the
maximum values of the components $u_{x}$ and $u_{y}$, although close, are not
strictly equal.

## References

  * Convery and Gadegaard (2019) N. Convery and N. Gadegaard, “30 years of microfluidics,” Micro and Nano Engineering 2, 76 (2019).
  * Whitesides (2006) G. M. Whitesides, “The origins and the future of microfluidics,” Nature 442, 368 (2006).
  * Beebe  _et al._ (2002) D. J. Beebe, G. A. Mensing, and G. M. Walker, “Physics and applications of microfluidics in biology,” Annu. Rev. Biomed. Eng 4, 261 (2002).
  * Stone  _et al._ (2004) H. A. Stone, A. D. Stroock, and A. Ajdari, “Engineering flows in small devices: Microfluidics Toward a Lab-on-a-Chip,” Annu. Rev. Fluid Mech. 36, 381 (2004).
  * Squires and Quake (2005) T. M. Squires and S. R. Quake, “Microfluidics: fluid physics at the nanoliter scale,” Rev. Mod. Phys. 77, 977 (2005).
  * Gu  _et al._ (2018) Y. Gu, V. Hegde, and K. J. M. Bishop, “Measurement and mitigation of free convection in microfluidic gradient generators,” Lab Chip 18, 3371 (2018).
  * Yoon  _et al._ (2005) S. K. Yoon, M. Mitchell, E. R. Choban, and P. J. A. Kenis, “Gravity-induced reorientation of the interface between two liquids of different densities flowing laminarly through a microchannel,” Lab Chip 5, 1259 (2005).
  * Selva  _et al._ (2012) B. Selva, L. Daubersies, and J.-B. Salmon, “Solutal convection in confined geometries: enhancement of colloidal transport,” Phys. Rev. Lett. 108, 198303 (2012).
  * Edwards  _et al._ (2018) 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. 121, 184501 (2018).
  * Li  _et al._ (2019) Y. Li, C. Diddens, P. Lv, H. Wijshoff, M. Versluis, and D. Lohse, “Gravitational Effect in Evaporating Binary Microdroplets,” Phys. Rev. Lett. 122, 114501 (2019).
  * Kang  _et al._ (2013) K. H. Kang, H. C. Lim, H. W. Lee, and S. J. Lee, “Evaporation-induced saline Rayleigh convection inside a colloidal droplet,” Phys. Fluids 25, 042001 (2013).
  * Pradhan and Panigrahi (2018) T. K. Pradhan and P. K. Panigrahi, “Convection inside condensing and evaporating droplets of aqueous solution,” Soft Matter 14, 4335 (2018).
  * Lee  _et al._ (2014) S. J. Lee, J. Hong, and Y.-S. Choi, “Evaporation-induced flows inside a confined droplet of diluted saline solution,” Langmuir 30, 7710 (2014).
  * Daubersies  _et al._ (2012) L. Daubersies, J. Leng, and J.-B. Salmon, “Confined drying of a complex fluid drop: phase diagram, activity, and mutual diffusion coefficient,” Soft Matter 8, 5923 (2012).
  * Loussert  _et al._ (2016) C. Loussert, A. Bouchaudy, and J.-B. Salmon, “Drying dynamics of a charged colloidal dispersion in a confined drop,” Phys. Rev. Fluids 1, 084201 (2016).
  * Inoue and Inasawa (2020) K. Inoue and S. Inasawa, “Drying-induced back flow of colloidal suspensions confined in thin unidirectional drying cells,” RSC Advances 10, 15763 (2020).
  * Pradhan and Panigrahi (2016) T. K. Pradhan and P. K. Panigrahi, “Evaporation-induced natural convection of a liquid slug of binary mixture inside a microchannel: effect of confinement,” Microfluidics and Nanofluidics 20, 115 (2016).
  * Savino and Monti (1996) R. Savino and R. Monti, “Buoyancy and surface-tension-driven convection in hanging-drop protein crystallizer,” J. Crystal Growth 165, 308 (1996).
  * Pradhan  _et al._ (2012) T. K. Pradhan, M. Asfer, and P. K. Panigrahi, “Droplet hydrodynamics during lysozyme protein crystallization,” Phys. Rev. E 86, 051602 (2012).
  * Apostolopoulou  _et al._ (2020) V. Apostolopoulou, N. Junius, R. P. Sear, and M. Budayova-Spano, “Mixing Salts and Poly(ethylene glycol) into Protein Solutions: The Effects of Diffusion across Semipermeable Membranes and of Convection,” Cryst. Growth Des. 20, 3927 (2020).
  * Pradhan and Panigrahi (2020) T. K. Pradhan and P. K. Panigrahi, “Suppressing internal convection of a droplet using confinement during protein crystallization,” J. Appl. Phys. 128, 084701 (2020).
  * Dunstan  _et al._ (2018) J. Dunstan, K. J. Lee, Y. Hwang, S. F. Park, and R. E. Goldstein, “Evaporation-driven convective flows in suspensions of nonmotile bacteria,” Phys. Rev. Fluids 3, 123102 (2018).
  * Hansen  _et al._ (2002) C. L. Hansen, E. Skordalakes, J. M. Berger, and S. R. Quake, “A robust and scalable microfluidic metering method that allows protein crystal growth by free interface diffusion,” Proc. Natl. Acad. Sci. USA 99, 16531 (2002).
  * Yamada  _et al._ (2016) 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 actionable microfilaments for biochemical assays, cell culture and organs-on-chip,” Lab Chip 16, 4691 (2016).
  * Hallez and Magnaudet (2008) Y. Hallez and J. Magnaudet, “Effects of channel geometry on buoyancy-driven mixing,” Phys. Fluids 20, 053306 (2008).
  * Shin  _et al._ (2004) J. O. Shin, S. B. Dalziel, and P. F. Linden, “Gravity currents produced by lock exchange,” J. Fluid Mech. 521, 1 (2004).
  * Martin  _et al._ (2011) J. Martin, N. Rakotomalala, L. Talon, and D. Salin, “Viscous lock-exchange in rectangular channels,” J. Fluid Mech. 673, 132 (2011).
  * Matson and Hogg (2012) G. P. Matson and A. J. Hogg, “Viscous exchange flows,” Phys. Fluids 24, 023102 (2012).
  * Séon  _et al._ (2007) T. Séon, J. Znaien, D. Salin, J. P. Hulin E. J., Hinch, and B. Perrin, “Transient buoyancy-driven front dynamics in nearly horizontal tubes,” Phys. Fluids 19, 123603 (2007).
  * Szulczewski and Juanes (2013) M. L. Szulczewski and R. Juanes, “The evolution of miscible gravity currents in horizontal porous layers,” J. Fluid Mech. 719, 82 (2013).
  * Tyn and Calus (1975) M. T. Tyn and W. F. Calus, “Temperature and Concentration Dependence of Mutual Diffusion Coefficients of Some Binary Liquid Systems,” J. Chem. Eng. Data 20, 310 (1975).
  * Crank (1975) J. Crank,  _The mathematics of diffusion_ (Oxford university press, 1975).
  * Xin and Le Quéré (2002) S. Xin and P. Le Quéré, “An extended Chebyshev pseudo-spectral benchmark for the 8:1 differentially heated cavity,” Int. J. Heat Mass Transfer 40, 981 (2002).
  * Xin  _et al._ (2008) S. Xin, J. Chergui, and P. Le Quéré, “3D spectral parallel multi-domain computing for natural convection flows,” in  _Parallel Computational Fluid Dynamics_ (2008) pp. 163–171.
  * Erdogan and Chatwin (1967) M. E. Erdogan and P. C. Chatwin, “The effects of curvature and buoyancy on the laminar dispersion of solute in a horizontal tube,” J. Fluid Mech. 29, 465 (1967).
  * Smith (1976) R. Smith, “Longitudinal dispersion of a buoyant contaminant in a shallow channel,” J. Fluid Mech. 78, 677 (1976).
  * Barton (1976) 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. 74, 81 (1976).
  * Godfrey (1980) J.S. Godfrey, “A numerical model of the James River estuary, Virginia, U.S.A.” Estuarine Coastal Mar. Sci. 11, 295 (1980).
  * Maclean and Alboussière (2001) D. J. Maclean and T. Alboussière, “Measurement of solute diffusivities. Part I. Analysis of coupled solute buoyancy-driven convection and mass transport,” Int. J. Heat Mass Transfer 44, 1639 (2001).
  * Salmon and Doumenc (2020) J.-B. Salmon and F. Doumenc, “Buoyancy-driven dispersion in confined drying of liquid binary mixtures,” Phys. Rev. Fluids 5, 024201 (2020).
  * Young and Jones (1991) W. R. Young and Scott Jones, “Shear dispersion,” Phys. Fluids A 3, 1087 (1991).
  * Chatwin and Sullivan (1982) P. C. Chatwin and P. J. Sullivan, “The effect of aspect ratio on longitudinal diffusivity in rectangular channels,” J. Fluid Mech. 120, 347 (1982).
  * Ajdari  _et al._ (2006) A. Ajdari, N. Bontoux, and H. A. Stone, “Hydrodynamic dispersion in shallow microchannels: the effect of cross-sectional shape,” Anal. Chem. 78, 387 (2006).
  * Shankar (1993) P. N. Shankar, “The eddy structure in Stokes flow in a cavity,” J. Fluid Mech. 250, 371 (1993).
  * Shankar  _et al._ (2002) P. N. Shankar, V. V. Meleshko, and E. I. Nikiforovich, “Slow mixed convection in rectangular containers,” J. Fluid Mech. 471, 203 (2002).
  * Joseph  _et al._ (1996) D. D. Joseph, A. Huang, and H. Hu, “Non-solenoidal velocity effects and Korteweg stresses in simple mixtures of incompressible liquids,” Physica D 97, 104 (1996).
  * Brenner (2005) H. Brenner, “Navier-Stokes revisited,” Physica A 349, 60 (2005).
  * Bruus (2007) H. Bruus,  _Theoretical Microfluidics_ , edited by Oxford Master Series in Physics (2007).