# Convective dissolution in confined porous media

Marco De Paoli1,2<EMAIL_ADDRESS>Christopher J. Howland1
<EMAIL_ADDRESS>Roberto Verzicco1,3,4<EMAIL_ADDRESS>Detlef
Lohse1<EMAIL_ADDRESS>1Physics of Fluids Group and Max Planck Center for
Complex Fluid Dynamics and J. M. Burgers Centre for Fluid Dynamics, University
of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands 2Institute of Fluid
Mechanics and Heat Transfer, TU Wien, 1060 Vienna, Austria 3Dipartimento di
Ingegneria Industriale, University of Rome “Tor Vergata”, 00133 Roma, Italy
4Gran Sasso Science Institute, 67100 L’Aquila, Italy 5Max Planck Institute for
Dynamics and Self-Organization, 37077 Gottingen, Germany

###### Abstract

We consider the process of convective dissolution in a homogeneous and
isotropic porous medium. The flow is unstable due to the presence of a solute
that induces a density difference responsible for driving the flow. The mixing
dynamics is thus driven by a Rayleigh-Taylor instability at the pore scale. We
investigate the flow at the scale of the pores using experimental
measurements, numerical simulations and physical models. Experiments and
simulations have been specifically designed to mimic the same flow conditions,
namely matching porosities, high Schmidt numbers, and linear dependency of
fluid density with solute concentration. In addition, the solid obstacles of
the medium are impermeable to fluid and solute. We characterise the evolution
of the flow via the mixing length, which quantifies the extension of the
mixing region and grows linearly in time. The flow structure, analysed via the
centre-line mean wavelength, is observed to grow in agreement with theoretical
predictions. Finally, we analyse the dissolution dynamics of the system,
quantified through the mean scalar dissipation, and three mixing regimes are
observed. Initially, the evolution is controlled by diffusion, which produces
solute mixing across the initial horizontal interface. Then, when the
interfacial diffusive layer is sufficiently thick, it becomes unstable,
forming finger-like structures and driving the system into a convection-
dominated phase. Finally, when the fingers have grown sufficiently to touch
the horizontal boundaries of the domain, the mixing reduces dramatically due
to the absence of fresh unmixed fluid. With the aid of simple physical models,
we explain the physics of the results obtained numerically and experimentally.
The solute evolution presents a self-similar behaviour, and it is controlled
by different length scales in each stage of the mixing process, namely the
length scale of diffusion, the pore size, and the domain height.

###### keywords:

convection, porous media, dissolution, mixing

##  1 Introduction

Natural convective flows take place when density differences within a fluid
domain drive the fluid motion. Natural convection is the dominant heat and
mass transport mechanism for many flows of practical interest, both in nature
(Hartmann et0al., 2001; Marshall & Schott, 1999; Cardin & Olson, 1994) and in
industrial applications (Bejan, 2013; Krepper et0al., 2002). One of the
particularly relevant ones in geophysics is the case of convection in porous
media, in which the fluid occupies the interstitial space of a solid porous
matrix. This configuration is typical of subsurface transport phenomena, and
is relevant for petroleum migration (Simmons et0al., 2001), water
contamination (LeBlanc, 1984; Van Der0Molen & Van0Ommen, 1988), underground
hydrogen storage (Zivar et0al., 2021; Krevor et0al., 2023), superficial
formations in dry salty lakes (Lasser et0al., 2021, 2023), and sea ice growth
(Wettlaufer et0al., 1997; Feltham et0al., 2006; Middleton et0al., 2016), to
name a few.

Recently, this problem has been subject of extended investigations due to the
implications it can bear for geological carbon sequestration (Huppert &
Neufeld, 2014; Emami-Meybodi et0al., 2015). In this process, carbon dioxide
(CO2) is injected into underground formations located one to three kilometres
beneath the earth surface, with the aim of permanent storage. After injection,
carbon dioxide remains in a supercritical state due to the high pressure at
the injections depths, but its density is lower compared to that of the
resident fluid (brine), and the injected volume of carbon dioxide migrates on
top of the brine layer (MacMinn et0al., 2012). This situation is undesired,
because it may cause a further migration of carbon dioxide to the upper
layers, eventually leading to a return of CO2 into the atmosphere. However,
carbon dioxide is partially soluble in brine, and when these fluids mix a
heavier solution (CO2 \+ brine) forms. When a sufficiently thick layer of this
denser mixture builds up, it may become unstable and form convective
instabilities (De0Paoli, 2021). These flow structures have a favourable effect
on the efficiency of the CO2 dissolution mechanism, since they contribute to
fluid mixing (i.e., CO2 permanently stored in brine) in a much more effective
manner compared to pure diffusion (Ennis-King et0al., 2005; Xu et0al., 2006).
On the other hand, the process is made complex by the presence of these non-
linear instabilities, and making reliable predictions of the dissolution
dynamics is a non-trivial task.

The dissolution of CO2 in brine is a problem characterised by an integral
length scale corresponding to the height of the reservoir, which may be of the
order of tens or hundreds of meters in geological formations of practical
interest (Huppert & Neufeld, 2014). Therefore, resolving all the flow scales,
from the pores to the reservoir scale, is not an option with the present
computational and experimental capabilities. A possible approach commonly
employed to revolve the large scales of the flow consists of analysing the
flow at the Darcy scale, i.e. the flow is not resolved within the pores, but
the quantities of interest (pressure, velocity, concentration, temperature)
are obtained as averaged over a representative volume containing many pores
(Nield & Bejan, 2017). The majority of the works on convective flows in porous
media currently available in literature refers to this case.

The mixing process of CO2 and brine is idealised considering an homogeneous
porous slab with a constant concentration difference applied at the horizontal
boundaries (Rayleigh-Bénard), or as a semi-infinite (one-sided) domain where
the concentration is fixed at top and no flux is applied at bottom boundary
(Hewitt, 2020). In the Rayleigh-Bénard-Darcy case, the system attains a
statistically steady-state, which has been predicted theoretically (Horton &
Rogers0Jr, 1945; Lapwood, 1948) and recently accurately described numerically
in two (Hewitt et0al., 2012; De0Paoli et0al., 2016; Wen et0al., 2018; Ulloa &
Letelier, 2022) and three dimensions (Hewitt et0al., 2014; De0Paoli et0al.,
2022b; Dhar et0al., 2022), also at high Rayleigh-Darcy numbers,
$\operatorname{\mathit{Ra}}^{*}$ (the Rayleigh-Darcy number indicates the
relative importance of convective and diffusive contributions in a flow
through a porous medium). The one-sided-Darcy configuration is characterised
by a transient behaviour (Riaz et0al., 2006; Fu et0al., 2013) consisting of
three main phases: (i) an initial diffusive regime in which the mixing layer
grows and becomes eventually unstable (De0Paoli et0al., 2017), (ii) a
convective phase characterised by a constant dissolution rate (Hidalgo et0al.,
2012), and (iii) a finite-size regime when the strength of convection reduces
gradually (Slim et0al., 2013).

The Rayleigh-Bénard-Darcy and the one-sided-Darcy systems have been well
characterised numerically and theoretically in case of thermal convection at
the Darcy scale, i.e. when the flow structures are large compared to the scale
of the pores (Hewitt, 2020). Numerical results suggest that a linear
relationship exists between the dimensionless mass transfer coefficient
(Sherwood number, $\operatorname{\mathit{Sh}}$) and the Rayleigh-Darcy number
($\operatorname{\mathit{Ra}}^{*}$), in both the Rayleigh-Bénard-Darcy case
(Hewitt et0al., 2012; Pirozzoli et0al., 2021) and during the convective regime
of the one-sided-Darcy system (Slim, 2014; De0Paoli et0al., 2017).
Experimental measurements in these configurations (Neufeld et0al., 2010;
Backhaus et0al., 2011; De0Paoli et0al., 2020; Brouzet et0al., 2022), however,
suggest a different qualitative and also quantitative behaviour compared to
the corresponding Darcy simulations, with a non-linear scaling of
$\operatorname{\mathit{Sh}}$ with $\operatorname{\mathit{Ra}}^{*}$. This
discrepancy is likely due to non-Darcy effects (Liang et0al., 2018), i.e. to
the pore-induced flow dynamics not captured by Darcy simulations. Therefore,
resolving the flow and the solute transport at the pore-scale is crucial to
make reliable models to incorporate in large-scale simulations and to predict
the underground fluids migration and mixing, and it represents the main
motivation for this work.

When a fluid flows through a matrix of solid obstacles, the fluid follows a
random-walk-type path (Woods, 2015). If the fluid is carrying a solute, in
addition to molecular diffusion, solute spreading may occur due to pore-scale
change of flow direction (mechanical dispersion), heterogeneities in the
aquifer (large-scale dispersion) or other mechanisms, such as dead-end pores
(anomalous dispersion). In this work, we will only refer to mechanical
dispersion, i.e. we assume that the medium is homogeneous and without dead-end
spaces. On the one hand, mechanical dispersion produces additional spreading
of solute and can be several orders of magnitude more effective than molecular
diffusion (Delgado, 2007). It has been also observed (Eckel et0al., 2022) that
the presence of the finger pattern and the counter-current flow structure
enhance the longitudinal spreading of the solute compared to unidirectional
dispersion of a single-solute plume. This additional spreading is responsible
for the reduction of the local density gradients, diminishing the strength of
convective motions and coupling convection and diffusion as mechanisms
controlling the mixing process. To quantify the relative importance of pore
structure, material properties and driving force on the overall heat or mass
transport process, pore-resolved convective flows have been recently employed
in the framework of thermal Rayleigh-Bénard convection, in both experiments
and simulations. Chakkingal et0al. (2019) observe that the flow structure and
the heat transfer coefficient are determined by the relative size of thermal
length scale (boundary layer thickness) and porous length scale (average pore
space). These properties determine the penetration of the plumes in the
boundary layer region, which is responsible for the heat or mass transfer
rate. In a complementary work, Ataei-Dadavi et0al. (2019) observed that while
at low Rayleigh numbers the transport mechanism is less efficient than in free
fluids Rayleigh-Bénard convection, at larger Rayleigh numbers the classical
$\operatorname{\mathit{Nu}}$ vs. $\operatorname{\mathit{Ra}}$ scaling
(Grossmann & Lohse, 2000, 2001) is recovered. The nature of this transition
has been investigated by Liu et0al. (2020). They used two-dimensional direct
numerical simulations to examine in detail the microscale flow field through a
bead pack. They observed that the transition between these two regimes is
controlled by two physical mechanisms induced by the porous matrix: (i) the
presence of obstacles makes the flow more coherent, with the correlation
between temperature fluctuation and vertical velocity enhanced and the
counter-gradient convective heat transfer suppressed, leading to heat transfer
enhancement; and (ii) the convection strength is reduced due the impedance of
the obstacle array, corresponding to heat transfer reduction. They observed
that the scaling crossover occurs when the thickness of the thermal boundary
layer is comparable to the averaged pore length scale. In addition to the
porous structure and the Rayleigh number, in case of thermal convection, the
boundary layer thickness and the heat transfer coefficient are determined also
by the value of thermal conductivity of the solid and liquid phases (Korba &
Li, 2022; Zhong et0al., 2023).

The discussed findings refer to systems in which the medium is permeable to
the scalar (heat or solute) transported. In the frame of solute dispersion,
two major differences arise compared to thermal convection, namely the solid
is usually impenetrable (over the flow time scales) to the solute, and the
ratio between momentum and solute diffusivity (Schmidt number,
$\operatorname{\mathit{Sc}}$) is typically about two orders of magnitude
larger than the ratio between momentum and heat diffusivity (Prandtl number,
$\operatorname{\mathit{Pr}}$). Gasow et0al. (2020, 2021) focused on solute
convection investigating solute transport at the pore-scale in Rayleigh-Bénard
configuration, with the aim of deriving corrections to be included in Darcy
models to account for the obstacles-induced solute dispersion. They observed
that the pore-induced dispersion, which may be as strong as buoyancy, affects
also the momentum transport and it is determined by two length scales (the
pore length scale and the domain height). Moreover, the dissolution
coefficient ($\operatorname{\mathit{Sh}}$) is observed to depend also on the
Schmidt number (Gasow et0al., 2022), in addition to the Rayleigh number and
the pore structure (Liu et0al., 2020). These numerical studies, which
represent a fundamental step to make large-scale predictions of dissolution in
porous media, are still computationally limited to two-dimensions, moderate
Rayleigh numbers and large porosity values. In comparison, experiments can
achieve larger Rayleigh numbers and values of porosity that are more
representative of geological formations. In contrast, in the instance of
solutal convection, the Rayleigh-Bénard and the semi-infinite configurations
are hard to tackle experimentally because of the limitations associated with
keeping the solute concentration constant at the boundaries. The porous
Rayleigh-Taylor instability, relatively less studied compared to the
corresponding Rayleigh-Bénard and one-sided configurations, is an excellent
candidate to overcome this obstacle, and it is the object of this work.

A porous Rayleigh-Taylor system consists of two miscible fluids of different
density initially arranged in an unstable configuration and immersed in a
porous matrix. After an initial diffusive phase (Wang et0al., 2016, 2018), the
flow is driven by convection, and due to the transient nature of this system,
the mixing rate of the two fluids varies in time. In geophysical applications,
it is crucial to determine the mixing rate of the involved fluids, to be able
to make reliable predictions on the evolution of the volume of subsurface
flow. The porous Rayleigh-Taylor system has been studied at the Darcy scale
with the aid of Hele-Shaw cells and Darcy simulations (De0Wit, 2020). The flow
behaviour is quantified by the mixing length, i.e. the vertical extension of
the mixing region. The growth rate of the mixing length is controlled by the
combined action of buoyancy (density difference) and drag (viscous dissipation
through the medium) (Boffetta & Musacchio, 2022; De0Paoli et0al., 2022a). In
the case of miscible fluids, the role of diffusion across the fluid-fluid
interface is also crucial (Gopalakrishnan et0al., 2017), as it weakens
convection by reducing local density gradients. As a result of this time-
dependent finger growth process, the mixing dynamics is also transient, and it
may be quantified via the mean scalar dissipation rate (Hidalgo et0al., 2012).
The mean scalar dissipation is particularly convenient because it can be
exactly related to other flow quantities, e.g. the averaged concentration
fluctuations, and it has been already described numerically at the Darcy scale
(De0Paoli et0al., 2019a). Also in the Rayleigh-Taylor case, however, the
behaviour of dispersion cannot be captured by Darcy simulations despite its
crucial role, as it can influence dramatically the onset of the gravitational
instabilities (Menand & Woods, 2005) and the mixing dynamics. In this flow
configuration, the full flow dynamics has been studied with the aid of three-
dimensional simulations by Sardina et0al. (2018), where the authors consider a
thermal convection at relatively large values of porosity ($0.6-1$). They
proposed a model to incorporate the effect of the medium within a friction
coefficient to be included in the Navier-Stokes equations. The case of solutal
convection in geological formations may be markedly different, since the
porosity is low (0.2 – 0.4) and the Schmidt numbers about two orders of
magnitude larger than in the thermal case, and we aim precisely at this gap.

We investigate a Rayleigh-Taylor instability in a saturated, homogeneous and
isotropic porous medium. We present original experiments and numerical
simulations, where we resolve the flow at the pore-scale at high Rayleigh-
Darcy numbers, high Schmidt numbers, and low values of porosity. Experiments
and simulations are specifically designed to reproduce the same fluid and
medium properties, namely linear density-concentration dependency, matching
porosities, and porous medium impermeable to the transported scalar (solutal
convection). In the experiments, porous media and fluid properties are varied,
and different flow regimes are observed, namely a Darcy-type flow (with the
flow structures larger than the pore length scale) and a diffusion-dispersion
regime, with the strength of mechanical dispersion being equivalent or
dominant with respect to molecular diffusion. Simulations are first validated
against the experimental measurements in terms of evolution of the mixing
length, and then used to quantify the evolution of the mixing rate, measured
by the mean scalar dissipation. Several dissolution regimes have been
identified: Initially, the flow is controlled by diffusion. This regime is
followed by a convection-dominated phase. The average concentration profiles
follow a self-similar behaviour, which we describe theoretically throughout
the flow dynamics. In the final regime finite-size effects reduce the solute
transport. Our experimental and numerical results are used to explain the
evolution of the flow with the aid of physically grounded models.

The paper is organised as follows. In §2 the problem is formulated. We
describe the experimental and the numerical setups in §3. We present our
results in terms of qualitative flow dynamics (§4), mixing length and
concentration profiles (§5), flow structure and wavenumber (§6). In §7 we
quantify the dissolution rate, which we describe by physical models in each of
the regimes identified. Finally, in §8 we summarise our findings and provide a
brief perspective on future research directions.

##  2 Problem formulation

The process of convective dissolution is studied here in the frame of the
Rayleigh-Taylor instability. It can be modelled as two layers of miscible
fluids having different density, initially separated by a flat interface, and
located in an accelerated reference frame (Boffetta & Mazzino, 2017). The
process is simulated in the context of porous media flows, mimicked with the
aid of a porous layer (height $H$, width $L$) made of spheres (diameter $d$).
The fluids fully saturate the medium, have same viscosity ($\mu$) but
different density, and are arranged in an unstable configuration, with the
heavy fluid (density $\rho_{0}$) lying on top of the lighter one (density
$\rho_{w}$). Therefore, the maximum density difference within the system is
$\Delta\rho=\rho_{0}-\rho_{w}$. The system, consisting of a porous slab
saturated with two fluids in an accelerated field (gravity acceleration, $g$),
is sketched in figure 1(a). The density difference is induced by the presence
of a solute, which is quantified by the solute concentration $C$, taking its
maximum at the upper layer ($C=C_{0}$) and minimum at the lower layer ($C=0$).
The reference frame ($x,z$) is defined as in figure 1(a) such that it is
centred at the mid height of the domain and $z$ is aligned with $g$ but in
opposite direction. We aim at mimicking a system with horizontal boundaries
($z=\pm H/2$) that are impermeable to both fluid and solute, and we consider

$\mathbf{u}\cdot\mathbf{n}=0\quad,\quad\frac{\partial C}{\partial n}=0\text{
,}$ (1)

with $\mathbf{n}$ the unit vector perpendicular to the boundary.

Figure 1:  Sketch illustrating the system considered. (a) Domain with explicit
indication of the boundary conditions (no flux of mass or solute through the
horizontal walls) and domain dimensions in horizontal ($L$) and vertical ($H$)
directions. The reference frame ($x,z$) as well as the initial position of the
interface (red dashed line) are indicated, with the heavy fluid (density
$\rho_{0}$, concentration $C=C_{0}$) initially lying on top of the lighter one
($\rho_{w}$, $C=0$). (b) In the experiments, a transparent medium consisting
of monodisperse beads and fluids of different colours are used. (c) In the
simulations, the geometry consists of an array of spheres (diameter $d$) fully
saturated with fluid. The fluid carrying the solute moves thought the spheres.

In the present flow configuration, the dimensionless parameters controlling
the system can be grouped in three-main categories: medium parameters (Darcy
number, porosity), fluid parameters (Schmidt number) and flow parameters
(Rayleigh, Rayleigh-Darcy, Peclet and Reynolds numbers).

We consider a simplified configuration in which the medium is homogeneous and
isotropic. Assuming the structure obtained from the sphere packing as an
isotropic and homogeneous medium, it can be fully described by two global
quantities, namely porosity and permeability. The porosity, $\phi$, represents
the ratio between the volume of fluid and the total volume (fluid + solid) of
the domain considered, and therefore it varies between $\phi=0$ (pure solid)
and $\phi=1$ (pure fluid). The permeability, $k$, quantifies the ability of
the porous matrix to allow a fluid to flow through it. For a given geometry of
the medium, the Darcy number

$\operatorname{\mathit{Da}}=\frac{k}{H^{2}}$ (2)

quantifies the relative dimension of the microscopic pore-scale ($\sqrt{k}$)
and the macroscopic length-scale ($H$) (Hewitt, 2020).

The dimensionless parameter that quantifies the fluid properties is the
Schmidt number, which is the ratio of momentum diffusivity (kinematic
viscosity, $\mu/\rho_{o}$) to mass diffusivity $D$,

$\operatorname{\mathit{Sc}}=\frac{\mu}{\rho_{0}D}\text{ .}$ (3)

The dimensionless flow parameters and the relevant flow scales are obtained by
combining domain, medium and fluid properties. A possible velocity scale of
the flow is the buoyancy velocity

$U=\frac{g\Delta\rho k}{\mu}\text{ ,}$ (4)

which is obtained at the equilibrium between driving forces ($gk\Delta\rho$)
and viscous dissipation through the medium ($\mu$).

Multiple length scales are effective in this problem. One can consider as a
reference length scale the distance $\ell$ over which advection and diffusion
balance (Slim, 2014)

$\quad\ell=\frac{\phi D}{U}\text{ .}$ (5)

Possible alternatives consists of the characteristic bead size (sphere
diameter, $d$) or the domain height ($H$). We will see that each of these
scales is relevant in different phases of the dissolution process. Solutal
convection in pure fluids is characterized by the competing effect of
convection (solute-induced density differences) and dissipation or diffusion,
respectively. The relative importance of these contributions is measured by
the concentration Rayleigh number based on the domain size
($\operatorname{\mathit{Ra}}$) or on the diameter of the spheres
($\operatorname{\mathit{Ra}}_{d}$),

$\operatorname{\mathit{Ra}}=\frac{g\Delta\rho H^{3}}{\mu
D}\quad\text{,}\quad\operatorname{\mathit{Ra}}_{d}=\frac{g\Delta\rho
d^{3}}{\mu D}=\operatorname{\mathit{Ra}}\left(\frac{d}{H}\right)^{3}\text{ ,}$
(6)

respectively. These parameters include convection and dissipation, but do not
consider the presence of the medium, which has a stabilizing effect on
convection due to the additional friction on the surface of the pores. The
ratio of the strength of these contributions is estimated by the Rayleigh-
Darcy number

$\operatorname{\mathit{Ra}}^{*}=\frac{g\Delta\rho kH}{\mu\phi
D}=\frac{UH}{\phi
D}=\frac{H}{\ell}=\frac{\operatorname{\mathit{Ra}}\operatorname{\mathit{Da}}}{\phi}.$
(7)

We remark that the concentration Rayleigh number (6) and the Rayleigh-Darcy
number (7) are linked to the porous medium properties via the Darcy number (2)
and the porosity. Finally, two more flow parameters are used to determine
whether the flow can be modelled as a Darcy flow or not. Following (Hewitt,
2020), the flow can be considered as a Darcy-type flow, if the length scale of
the flow structures is much larger than the representative volume over which
the quantities are averaged. It is obtained for (i) viscous forces dominating
over inertia at the pore-scale, and (ii) length scale of the convective flow
large compared to the pore size. These conditions are fulfilled if

$\Rey=\frac{\operatorname{\mathit{Ra}}^{*}\operatorname{\mathit{Da}}^{1/2}}{Sc}\ll
1\quad,\quad\operatorname{\mathit{Pe}}=\operatorname{\mathit{Ra}}^{*}\operatorname{\mathit{Da}}^{1/2}\ll
1,$ (8)

with $\Rey$ and $\operatorname{\mathit{Pe}}$ the pore-scale Reynolds number
and the Peclet number, respectively. In this study, only a few experiments
(and no simulations) fall in the Darcy case, and the relative flow dynamics
will be discussed later in §4. Note that in this definition of
$\operatorname{\mathit{Pe}}$ it is assumed that the pore-scale length used as
a length-scale for $\operatorname{\mathit{Pe}}$ is $\sqrt{k}$. An alternative
choice consists of using $d$, which would produce larger values of
$\operatorname{\mathit{Pe}}$ (by a factor of approx. 7.5 in this
configuration).

The flow scales of the experiments are listed in table 3.1, while the
dimensionless parameters corresponding to present experiments and simulations
is reported in table 3.2.

##  3 Methodology

We investigate convective mixing in confined porous media. The flow driving
force consists of the density differences induced by the presence of a solute.
We consider two miscible fluids characterized by a linear dependency of
density with concentration. The fluids are immersed in a fully saturated,
homogeneous and isotropic porous medium. Laboratory experiments and numerical
simulations are used to investigate the problem. In §3.1 we present the
experimental setup, consisting of a bead pack and an optical measurement
system. Pore-resolved two-dimensional simulations, where circular impermeable
obstacles are employed to mimic the solid matrix of the porous medium, are
discusses in §3.2.

###  3.1 Laboratory Experiments

The laboratory experiments are performed with the aid of a thick Hele-Shaw
cell filled with monodisperse beads and saturated with two fluids of different
densities in an unstable configuration. The parameters that can be varied are
the density difference $\Delta\rho$ between the fluids, and the diameter $d$
of the beads. Combining these parameters one can determine the flow reference
scales, namely $\ell$ and $U$. The experiments performed are listed in table
3.1.

Figure 2:  (a) Schematic representation of the Hele-Shaw cell. The cell is
filled with beads and different fluids from the upper and lower valves. The
gate keeps the fluids initially separated. The reference frame of the lab
($x,y,z$) is shown. The area of the gate (highlighted in green) and the
measurement region (highlighted in red) are discussed in detail in the panels
(b) and (c). (b) Side view of the gate mechanism. The gate is removed and the
fluids can mix. A system of seals is used to prevent mixing of the fluids
before the gate is open. After the gate is extracted, the cell is rotated
about the $x$ axis (i.e., upside down) as indicated by the red arrow. (c)
Front view of the measurement region. After the cell is rotated, the fluids
are in an unstable configuration (heavy solution on top of the lighter fluid)
and the recording starts.

In the following, we will discuss the experimental setup, the measurement
procedure, the fluid properties and the porous medium properties. The employed
Hele-Shaw cell used is sketched in figure 2(a). It consists of a hexagonal
container with uniform thickness. The shape is designed to simplify the fluid
extraction/injection phases. The hexagonal walls (made of PMMA layers,
thickness 8 mm) are transparent to allow optical access to the flow, and are
kept in place by a thick frame (25 mm) and bolts. A gasket (2 mm thick) is
placed between the frame and the walls to prevent fluid leakage. The cell is
initially filled with monodisperse glass spheres, then it is vibrated to
consolidate the medium, and finally the two fluids of different density are
injected via the upper and lower valves. Note that a gate (thickness 1 mm,
material polyethylene terephthalate - PET), as indicated in the side view in
figure 2b), initially divides the upper and the lower sides of the cells to
prevent fluids from mixing. A 1 mm high housing, equipped with a rubber seal
all along its length, is obtained on the front wall, and the gate fits in. The
gate can slide through an opening located on the back side of the cell, which
has additional seals hold in place by a plastic frame.

When the medium is consolidated and the fluids are injected, the cell is
placed in a stable configuration (with light fluid on top of the heavy one).
The gate is later removed and the cell rotates about the $x$ axis (red arrow
in figure 2b) to turn in the initial unstable condition (heavy fluid on top of
the light one) chosen to initialize the experiment. The measurements are
performed only in the central portion of the domain, indicated in red in
figure 2(a) and discussed in figure 2(c), which has size $H\times L\times
B=200\times 200\times 28.5$ mm3. The cell is illuminated on one side by a LED
light (Phlox LEDW-BL $300\times 300$ mm2) and on the opposite side a high-
resolution camera (Nikon D850 with lenses AF Nikkor, 50 mm, 1:1.8D) records
the evolution of the flow at an acquisition rate that varies between 0.05 and
2 frames per second.

To allow here a reliable comparison of experimental results against the direct
numerical simulations, the fluids employed in the experiments have been
specifically selected because of the linear dependency of the density of the
solution with respect to the solute concentration (Slim et0al., 2013; De0Paoli
et0al., 2022a). This condition is indeed met not only in the simulations
presented here, but also in most of the numerical works available in the
literature.

The employed fluids are water and an aqueous solution of KMnO4 (Potassium
Permanganate, Thermo Scientific, ACS reagent). We consider that the dynamic
viscosity, $\mu=9.2\times 10^{-4}$ Pa$\cdot$s, is constant and independent of
the solute concentration (Slim et0al., 2013). Similarly, we assume that the
diffusion coefficient is not sensibly affected by solute concentration, and
corresponds to $D=1.65\times 10^{-9}$ m2/s. While water density ($\rho_{w}$)
is nearly constant among the experiments (it is only dependent on temperature,
$\vartheta$), the density of the aqueous solution of KMnO4 can be varied by
changing the solute concentration in the aqueous solution, $C$, which is used
to control the density difference between the heavy and the light fluid. The
mass fraction of the solution is defined as

$\omega(C)=\frac{C}{\rho(C)}\text{ .}$ (9)

The respective dependency of density, mass fraction and concentration can
easily be determined. The density of the mixture, $\rho(C)$, can be well
approximated by a linear function of the solute concentration, i.e., it meets
the desired condition:

$\rho=\rho_{0}\biggl{[}1+\frac{\Delta\rho}{\rho_{0}C_{0}}\bigl{(}C-C_{0}\bigr{)}\biggr{]}\textit{
.}$ (10)

The concentration-density profiles as well as additional details are reported
in Appendix A.1.

With the aim of mimicking a homogeneous and isotropic porous medium, we fill
the cell with monodisperse spheres having diameter $d$, with 1 mm $\leq d\leq
4$ mm. Provided that the spheres are monodisperse, the diameter of the beads
and the porosity of the medium are the two parameters that determine the
medium property, i.e., the permeability. In the following, we will discuss
bead size, medium porosity and permeability. Again, a summary of all the
experimental parameters considered is reported in table 3.1.

The porosity of the medium indicates the ratio between the volume of fluid
used to fill the cell and the total cell volume (fluid and beads). We measure
the porosity by comparing the volume of water required to fill the cell with
and without the beads. The preparation of the medium is crucial in determining
the cell porosity and permeability. In this work, the cell is vibrated before
injecting the fluid so that the medium consolidates. Following this procedure,
the beads form a close random packing and the expected value of porosity in
case of monodisperse spheres is $\phi=0.359-0.375$ (Dullien, 1991; Haughey &
Beveridge, 1969). The values of porosity measured experimentally are
$\phi=0.37$ for all nominal diameters considered, except $d=3.00$ mm, in which
the value of porosity measured is slightly lower ($\phi=0.35$). This
difference can be possibly attributed to the lower quality of the beads with
$d=3.00$ mm, which have a wide distribution of diameters (see Appendix A.2).
Indeed, the more dispersed the diameters, the lower the value of porosity that
can be achieved.

Name |  | Fluid |  | Medium |  | Flow
---|---|---|---|---|---|---
|  | $\vartheta$ | $\omega_{M}$ | $\Delta\rho$ |  | $d$ | $\phi$ | $k$ |  | $\ell$ | $U$
|  |  [∘C] | [-] |  [kg/m3] |  | [mm] | [-] |  [m2] |  | [m] | [m/s]
E1 |  | 19.8 | 0.00112 | 0.877 |  | 1.00 | 0.37 |  7.09$\times 10^{-10}$ |  |  9.205$\times 10^{-5}$ |  6.632$\times 10^{-6}$
E2 |  | 19.6 | 0.00225 | 1.760 |  | 1.00 | 0.37 |  7.09$\times 10^{-10}$ |  |  4.588$\times 10^{-5}$ |  1.331$\times 10^{-5}$
E3 |  | 20.0 | 0.00450 | 3.528 |  | 1.00 | 0.37 |  7.09$\times 10^{-10}$ |  |  2.289$\times 10^{-5}$ |  2.667$\times 10^{-5}$
E4 |  | 20.0 | 0.00895 | 7.035 |  | 1.00 | 0.37 |  7.09$\times 10^{-10}$ |  |  1.148$\times 10^{-5}$ |  5.319$\times 10^{-5}$
E5 |  | 20.0 | 0.00116 | 0.903 |  | 1.75 | 0.37 |  2.17$\times 10^{-9}$ |  |  2.921$\times 10^{-5}$ |  2.090$\times 10^{-5}$
E6 |  | 19.6 | 0.00225 | 1.760 |  | 1.75 | 0.37 |  2.17$\times 10^{-9}$ |  |  1.498$\times 10^{-5}$ |  4.075$\times 10^{-5}$
E7 |  | 19.7 | 0.00450 | 3.521 |  | 1.75 | 0.37 |  2.17$\times 10^{-9}$ |  |  7.489$\times 10^{-6}$ |  8.152$\times 10^{-5}$
E8 |  | 19.5 | 0.00894 | 7.014 |  | 1.75 | 0.37 |  2.17$\times 10^{-9}$ |  |  3.759$\times 10^{-6}$ |  1.624$\times 10^{-4}$
E9 |  | 20.1 | 0.00111 | 0.868 |  | 3.00 | 0.35 |  5.07$\times 10^{-9}$ |  |  1.229$\times 10^{-5}$ |  4.699$\times 10^{-5}$
E10 |  | 19.6 | 0.00235 | 1.837 |  | 3.00 | 0.35 |  5.07$\times 10^{-9}$ |  |  5.812$\times 10^{-6}$ |  9.936$\times 10^{-5}$
E11 |  | 20.0 | 0.00453 | 3.547 |  | 3.00 | 0.35 |  5.07$\times 10^{-9}$ |  |  3.009$\times 10^{-6}$ |  1.919$\times 10^{-4}$
E12 |  | 19.2 | 0.00905 | 7.098 |  | 3.00 | 0.35 |  5.07$\times 10^{-9}$ |  |  1.504$\times 10^{-6}$ |  3.840$\times 10^{-4}$
E13 |  | 19.7 | 0.00112 | 0.872 |  | 4.00 | 0.37 |  1.13$\times 10^{-8}$ |  |  5.791$\times 10^{-6}$ |  1.054$\times 10^{-4}$
E14 |  | 19.7 | 0.00225 | 1.760 |  | 4.00 | 0.37 |  1.13$\times 10^{-8}$ |  |  2.867$\times 10^{-6}$ |  2.129$\times 10^{-4}$
E15 |  | 19.9 | 0.00451 | 3.529 |  | 4.00 | 0.37 |  1.13$\times 10^{-8}$ |  |  1.430$\times 10^{-6}$ |  4.269$\times 10^{-4}$
E16 |  | 19.9 | 0.00892 | 7.007 |  | 4.00 | 0.37 |  1.13$\times 10^{-8}$ |  |  7.203$\times 10^{-7}$ |  8.475$\times 10^{-4}$
Table 1:  Dimensional parameters of the performed experiments. The domain size
is constant ($L=H=200$0mm, $B=28.5$0mm). Fluid temperature ($\vartheta$),
solute mass fraction ($\omega_{M}$), and density contrast ($\Delta\rho$) are
reported. Fluid viscosity and diffusivity are $\mu=9.2\times 10^{-4}$
Pa$\cdot$s and $D=1.65\times 10^{-9}$0m2/s, respectively. Water density is
assumed $\rho_{0}=10^{3}$0kg/m3. Bead diameter ($d$), medium porosity
($\phi$), and permeability ($k$) are indicated. Flow length scale $\ell$ and
velocity scale $U$, defined as in (5) and (4) respectively, are also reported.

The permeability is inferred from the Kozeny-Carman correlation, i.e.

$k=\frac{d^{2}}{36k_{C}}\frac{\phi^{3}}{(1-\phi)^{2}}$ (11)

where $k_{C}$ is the Carman constant. This phenomenological correlation is
obtained for creeping flow, and it is found to be independent of the particle
shape (Dullien, 1991) (for non-spherical particles, $d$ is the equivalent
diameter). The Carman constant originally proposed within the Kozeny-Carman
formulaton (monodisperse spheres) is $k_{C}=4.8\pm 0.3$, usually approximated
to 5, which gives $36k_{C}=180$. At a later time, Ergun proposed the Blake-
Kozeny formulation, in which the coefficient is $36k_{C}=150$. Both
formulations are based on fitting of experimental data for flow through
granular beds. Other formulations have been proposed to improve the
predictions at low or high values of porosity, as well as to include the
effect of the Reynolds number. A detailed review on the possible values of
$k_{C}$ is provided by Xu & Yu (2008), where a more general phenomenological
formulation is also proposed. For the specific case of monodisperse spheres
randomly packed, Zaman & Jalali (2010) have shown that the Kozeny-Carman
formulation0(11) with $k_{C}=5$ provides good results, and this is the
correlation we employ.

The flow is recorded by the camera. The raw images of the measurement region
are analysed to quantify the flow evolution. An example is shown in figure03,
where an instantaneous intensity field obtained from experiment E12 is
discussed. The measurement region (see figure02c) consists of the central
squared portion of the cell, having size $H\times H=(200$0mm$)^{2}$. The
denser solution, initially on top, is much darker than the lighter one
(transparent). The thickness of the cell is large, and the light intensity
detected at the upper layer is weakly evolving during the experiment.
Therefore, for better accuracy and due to symmetry of the system, we consider
only the lower portion of the cell for the experimental measurements ($z/H\leq
0$).

Figure 3: Processing of images. (a)0Normalised intensity profile $I_{1}$,
defined in equation0(12), shown over the entire measurement region for the
experiment E12. The initial position of the interface ($z/H=0$, red solid
line) is reported. The instantaneous position of the interface, defined as an
iso-contour of $I_{1}$, is also shown (blue line). (b)0Corrected intensity
profile $I_{2}$, defined as in equation0(13). The corrected intensity $I_{2}$
is lower than 1 within the mixing region. (c)0Horizontally-averaged intensity
profiles. The lower edge of the mixing region, determined by the position
where $\overline{I}_{2}$ is lower than a given threshold value, is reported
(dashed red line). This is related to the mixing length $h$, discussed in
detail in §5.

The pre-processed fields are analyzed to produce intensity profiles and infer
the time-dependent evolution of the interface. The pre-processing consists of
several steps, illustrated in figure03. The initial light intensity
distribution $I(x,z,t=0)$, obtained from the raw images, is used to compute
and store the initial mean light intensity values of the high ($I_{H}$,
$z/H>0$) and low ($I_{L}$, $z/H<0$) fluid density layers. The normalised light
intensity field

$I_{1}(x,z,t)=\frac{I(x,z,t)-I_{H}}{I_{L}-I_{H}}.$ (12)

is suitable to visualise the instantaneous flow configuration. It is used, for
instance, to identify the instantaneous position of the interface of the
mixing region (figure03a). However, we observe in figure03(c) that the
horizontally-averaged intensity profile $\overline{I}_{1}$ varies smoothly
within the mixing region, and therefore it is not a good indicator to
determine the edge of the interface. We introduce the corrected intensity
(figure03b), $I_{2}$, defined as

$I_{2}(x,z,t)=2\left[\frac{1}{2}+I_{1}(x,z,t)-\widehat{I}_{1}(x,z,0)\right],$
(13)

with $\widehat{I}_{1}(x,z,0)$ the initial intensity field $I_{1}(x,z,t=0)$
computed with a moving average (squared window of size 100pixel). We observe
that $I_{2}$ is lower than 1 only within the mixing region and this property,
which is also clear from the horizontally-averaged profiles in figure03(c),
makes $I_{2}$ a more reliable observable to quantify the extension of the
mixing region (red dashed line).

###  3.2 Numerical simulations

Figure 4:  Schematic of the porous medium layout in the simulations. $(a)$
Regular hexagonal distribution of solid circles used in the simulations. $(b)$
An example of random perturbations to the regular pattern. Grey areas mark the
possible area in which each solid circle can move without overlapping another.
$(c)$ Schematic of the total domain for $H/d=70$, with solid particles in
grey, and the fluid domain coloured according to the initial concentration.

In the numerical part of this study, we solve the Navier–Stokes equations for
momentum, subject to the Oberbeck–Boussinesq approximation. This means that
variations in the density are only significant in the buoyancy term. We assume
the flow to be incompressible and impose
$\boldsymbol{\nabla}\cdot\boldsymbol{u}=0$ on the velocity field
$\boldsymbol{u}$. We consider variations in density to be linearly related to
the concentration field $C$, which itself satisfies an advection-diffusion
equation. We therefore consider the following governing equations

$\displaystyle\partial_{t}\boldsymbol{u}+(\boldsymbol{u}\cdot\boldsymbol{\nabla})\boldsymbol{u}$
$\displaystyle=-\rho_{0}^{-1}\boldsymbol{\nabla}p+\nu\nabla^{2}\boldsymbol{u}-g\beta
C\hat{\boldsymbol{z}},$ (14)
$\displaystyle\partial_{t}C+(\boldsymbol{u}\cdot\boldsymbol{\nabla})C$
$\displaystyle=D\nabla^{2}C,$ (15)

where $\rho_{0}$ is a reference density, $\nu=\mu/\rho_{0}$ is the kinematic
viscosity, $D$ the solutal diffusivity, $g$ gravitational acceleration, and
$\beta$ the solutal contraction coefficient describing the linear relationship
between density and concentration. The pressure $p$ satisfies a Poisson
equation such that a divergence-free velocity field is ensured.

We solve the equations (14)-(15) in a two-dimensional domain of height $H$ and
of width $\sqrt{3}H$. Periodic boundary conditions are imposed in the
horizontal ($x$) direction for all flow variables. At the top and bottom walls
($z=\pm H/2$), we impose the boundary conditions0(1), i.e. zero mass flux of
solute ($\partial_{z}C=0$) and no penetration ($w=0$, being $w$ the vertical
velocity). In addition, since the pore-scale flow is resolved, also the no-
slip condition applies along these walls ($\boldsymbol{u}=\boldsymbol{0}$). In
the following subsections, we describe the properties of the solid porous
matrix that occupies a portion of the simulation domain, as well as details of
the numerical implementation for the flow solver and the conditions at the
liquid-solid boundaries.

The numerical simulations are designed to match the porosity of the
experiments, namely $\phi=0.37$. To consider a two-dimensional setup most
similar to the monodisperse spherical bead pack of the experiments, we
construct the solid phase in the simulations from circles of a given diameter
$d$. Most past studies of a similar configuration (e.g. Sardina et0al., 2018;
Liu et0al., 2020) that explicitly resolve the pore-scale dynamics with a
liquid-solid boundary are performed at a higher porosity, allowing for a wide
range of configurations for the solid phase. Since we aim to match a low
experimental porosity of $\phi=0.37$, we prescribe a hexagonal arrangement of
the circular beads, as shown in figure 4(a), which allows for free percolation
of the fluid in 2D at these low porosities. Such perfectly regular
arrangements are not representative of the porous matrix in the experiments,
so we also repeat our simulations in domains in which random shifts from this
hexagonal arrangement are made to the positions of the solid circles. An
example of these random shifts is shown in figure 4(b). To prevent regions of
trapped fluid, we limit the random perturbations to the grey areas in that
schematic, such that the (black) solid regions do not overlap.

As discussed above in §3.1, the permeability $k$ is key to understanding the
effect of the porous medium properties on the flow. For example, the velocity
scale $U$ of (4) relies on the balance between buoyancy, permeability, and
viscosity. While the determination of the permeability for three-dimensional
arrays of spheres is well studied, for two-dimensional flows (array of
cylinders with infinite length) the situation has been investigated less. By
definition, the value of permeability would be determined by measuring the
pressure drop across the medium for different flow rates. Happel & Brenner
(2012) suggest that $k_{C}=5$ holds also for two-dimensional media. They
considered a flow perpendicular to an array of cylinders (indicated as
perpendicular flow in Xu & Yu, 2008) and observed that for $0.25<\phi<0.55$,
the Carman constant can very well approximated as $k_{C}=5$. Therefore, also
in the two-dimensional case, we will assume that equation0(11) with $k_{C}=5$
applies.

We use the highly-parallelised AFiD (Advanced Finite-Difference) code to
perform our simulations. The governing equations (14)-(15) are solved using
second-order central finite differences to compute spatial derivatives, with
time-stepping performed using an implicit Crank-Nicolson method for the
vertical diffusive terms ($\partial_{zz}$) and a third-order explicit Runge-
Kutta scheme for all other terms. A fractional-step method is used to impose
the divergence-free condition at each time step, where a Poisson equation is
solved for the pressure. Further details of the numerical scheme can be found
in Verzicco & Orlandi (1996) and van der Poel et0al. (2015). A multiple-
resolution method is applied to the concentration field for accurate
simulation at low diffusivities, following Ostilla-Monico et0al. (2015). Cubic
Hermite interpolation is used to interpolate the concentration field to the
velocity grid for the buoyancy forcing, and to interpolate the velocity field
to the refined grid for advection of the concentration field.

The solid phase is handled using the immersed boundary method (Verzicco,
2023). We follow the direct forcing approach of Fadlun et0al. (2000) such that
zero velocity is imposed in the solid during the implicit step of the
numerical solution. In this approach, linear interpolation is used at the
first grid nodes in the fluid from the solid boundary, so that the interface
is captured more accurately than the grid resolution. Similarly for the
concentration field, we ensure zero normal gradient at the immersed boundary
by specially treating these boundary nodes.

The collection of performed simulations are listed in table03.2. For each set
of parameters, we perform three simulations: one with the regular hexagonal
packing defining the centre positions of the solid circles, and two with
(distinct) random perturbations to this arrangement. Each simulation is
initialised with zero velocity, and an initial concentration field of

$\left.C(z)\right|_{t=0}=C_{0}\left[\mathscr{H}(z)+A(z)W\right],\quad
A(z)=0.1\mathrm{sech}^{2}\left(100\,z/H\right)$ (16)

where $\mathscr{H}$ is the Heaviside function, and $W$ is a white noise random
variable taking values between $-1$ and $1$, producing a region of random
perturbations of width $H/100$ across the interface. Uniform grid spacing is
used in all directions, with a resolution of at least 32 grid points per solid
diameter for the velocity grid and a resolution 4 times that for the refined
concentration grid. The largest computational grids are used for simulations
S10-S12, where $H/d=70$, with the base grid at a resolution of $4096\times
2365$ and the refined grid at a resolution of $16384\times 9459$. This
resolution is noticeably higher than in previous studies (e.g. Sardina et0al.,
2018), and allows us to accurately simulate the small-scale structures arising
from buoyancy-driven flows at high $Sc$.

Table 2: Dimensionless parameters of the performed experiments (E#) and the
simulations (S#). The domain aspect ratio of the experiments and the
simulations is 1 and $\sqrt{3}$, respectively.

S1 & 17 0.37 100 5.268$\times 10^{8}$ 1.000$\times 10^{5}$ 3.334$\times
10^{3}$ 5.102 0.0510 S2 17 0.37 100 1.666$\times 10^{9}$ 3.162$\times 10^{5}$
1.054$\times 10^{4}$ 16.135 0.1614 S3 17 0.37 100 5.268$\times 10^{9}$
1.000$\times 10^{6}$ 3.334$\times 10^{4}$ 51.024 0.5102 S4 35 0.37 100
4.214$\times 10^{9}$ 1.000$\times 10^{5}$ 6.669$\times 10^{3}$ 5.102 0.0510 S5
35 0.37 100 1.333$\times 10^{10}$ 3.162$\times 10^{5}$ 2.109$\times 10^{4}$
16.135 0.1614 S6 35 0.37 100 4.214$\times 10^{10}$ 1.000$\times 10^{6}$
6.669$\times 10^{4}$ 51.024 0.5102 S7 52 0.37 100 1.422$\times 10^{10}$
1.000$\times 10^{5}$ 1.000$\times 10^{4}$ 5.102 0.0510 S8 52 0.37 100
4.498$\times 10^{10}$ 3.162$\times 10^{5}$ 3.163$\times 10^{4}$ 16.135 0.1614
S9 52 0.37 100 1.422$\times 10^{11}$ 1.000$\times 10^{6}$ 1.000$\times 10^{5}$
51.024 0.5102 S10 70 0.37 100 3.372$\times 10^{10}$ 1.000$\times 10^{5}$
1.334$\times 10^{4}$ 5.102 0.0510 S11 70 0.37 100 1.066$\times 10^{11}$
3.162$\times 10^{5}$ 4.218$\times 10^{4}$ 16.135 0.1614 S12 70 0.37 100
3.372$\times 10^{11}$ 1.000$\times 10^{6}$ 1.334$\times 10^{5}$ 51.024 0.5102

##  4 Flow dynamics

We now present the results of the experiments and simulations. The experiments
are performed by changing the bead diameter ($d$) and the density difference
($\Delta\rho$). Simulations are performed for different values of diameter-
based Rayleigh number $(\operatorname{\mathit{Ra}}_{d})$ and domain to bead
size $(H/d)$.

Under certain flow conditions, a fluid flow through a porous medium may be
considered as a Darcy flow. In a Darcy flow, the length scale of the flow
structures is much greater than the representative volume over which the
quantities are averaged (Hewitt, 2020). This representative volume typically
includes a number of solid particles and the interstitial fluid. Darcy
conditions are met when (i)0the flow is controlled by viscous forces at the
pore scale $(\Rey\ll 1)$, and (ii)0the length scale of the convective flow is
large compared to the pore size $(\operatorname{\mathit{Pe}}\ll 1)$. One can
observe from table03.2 that only few experiments (E1-E3, E5) fall in the Darcy
case. A qualitative observation of this result is provided by looking at the
raw images in figure05 and 06.

Figure 5:  Flow evolution (raw images of laboratory experiments) for different
bead size, $d$, and same density difference, $\Delta\rho\approx 7$0kg/m3. The
permeability is increased while the driving force remains unchanged. The flow
structure changes remarkably from (a) to (d), due to the change in the size of
the fingers relative to the beads. For visualisation purposes, only the
central lower portion of the domain is shown, approximately corresponding to
$1/8\leq x/H\leq 7/8$ and $z/H\leq 0$ (length of scale bar is 2 cm).

In figure05, we consider the variation of the flow topology for the same
driving force ($\Delta\rho\approx 7$0kg/m3) and different values of
permeability (i.e., different $d$). Both $\Rey$ and
$\operatorname{\mathit{Pe}}$ are sensitive to variations of the diameter of
the beads. However, $\Rey$ remains reasonably lower than unity for E4, E8, E12
and E16, whereas the Peclet number changes remarkably up to $O(10^{2})$. This
reflects the fact that the length scale of the convective flow is of the same
order or smaller than the pore size. The transition is apparent when going
from E4 in figure05(a), where the fingers width is at least a few diameters,
to E12 in figure05(c), where one single plume penetrates and branches into the
interstitial space. In E16, shown in figure05(d), the permeability is
considerably larger than in previous cases. As a result, the driving force is
extremely vigorous compared to the viscous drag, with the velocity $U$
approximately twice larger than in E12 (see table03.1). With
$\operatorname{\mathit{Pe}}\approx 150$, the solute spreads quickly also in
the direction perpendicular to the view, which makes the light intensity field
blurry.

Figure 6:  Flow evolution (raw images of laboratory experiments) for the same
bead size, $d=1.75$0mm, and different density differences, $\Delta\rho$. The
driving force is increased while the permeability remains unchanged. The flow
structure changes remarkably from (a) to (d), due to the change in the size of
the fingers relative to the beads. For visualisation purposes, only the
central lower portion of the domain is shown, approximately corresponding to
$1/8\leq x/H\leq 7/8$ and $z/H\leq 0$ (length of scale bar is 2 cm).

In figure06, we consider the opposite scenario, i.e., for the same
permeability ($d=1.75$0mm) we vary the driving force (different $\Delta\rho$).
In this case, there is an increase in $\Rey$, which however remains
considerably lower than 1. The Peclet number increases of about one order of
magnitude between E5 and E8, and we can appreciate this smooth transition from
the intensity fields. The $\operatorname{\mathit{Pe}}$ of E5 is just above
unity and the width of the flow structures, shown in figure06(a), is about
$5d-10d$. Increasing the density difference makes the structures progressively
smaller (E6, E7), and eventually the fingers propagate as the thin filamentary
plumes (E8) shown in figure06(d). Although the characteristic size of these
plumes is comparable to the pore diameter, splitting in multiple branches is
still observed.

Figure 7:  Snapshots of the concentration field from simulations S10, S12, and
S6. In the first two simulations, $H/d=70$ is fixed, and only the Rayleigh
number $Ra_{d}$ is varied from $10^{5}$ to $10^{6}$ (left to centre). In the
third simulation (S6, right column), $Ra_{d}=10^{6}$ matches the centre column
but $H/d$ is decreased to 35, so the obstacles are larger compared to the
domain height. Each snapshot is presented at a given fraction of the timescale
$H/U$, with time increasing top to bottom (as specified in the left column).

We have shown that the experimental results may present in some cases features
that can be captured by a Darcy model when $Pe$ is low. By contrast, the
Peclet number for each numerical simulation satisfies $Pe\geq 5$, producing a
flow that does not fulfil the Darcy assumptions outlined in Eq.0(8). Indeed,
thin fingers penetrating the pore-space are always observed in the
simulations, as shown in figure 7 through a series of snapshots of the
concentration field for various pore-scale Rayleigh numbers $Ra_{d}$. In each
case, the early-time snapshots highlight an initial instability on the scale
of the beads, with thin plumes growing from the interface at $z=0$. As time
goes on, larger structures begin to develop in the mixing layer, with coherent
fingers spanning the centre of the domain. Even as these large-scale fingers
grow, the dynamics at the tips of the fingers continues to consist of thin
percolating structures, similar to what is observed in the experiments in
figures 6(b)-(d).

Comparing the first two columns in figure07 with each other, the main effect
of changing $Ra_{d}$ is in the lateral diffusion of the concentration field.
The late-time concentration observed in panel (j) appears somewhat smeared
out, with smoother gradients compared to the equivalent panel at higher
$Ra_{d}$, panel (k). The same Rayleigh number but a different geometry of the
medium are considered in the central and right columns. The minimum structure
size is set by the pore-space, and the width of the fingers (relative to $H$)
increases from $H/d=70$ (central column) to the $H/d=35$ (right column). In
all the simulations, the concentration field is strongly influenced by the
hexagonal lattice structure of the medium, with fingers percolating aligned at
an angle. Although figure 7 only features simulations with regular bead
patterns, we shall show in the following sections that quantitative measures
of the mixing are not significantly affected by the random obstacle
perturbations shown in figure 4. Furthermore, when time is scaled with $H/U$
as in figure 7, the size of the mixing layer appears consistent across all
three simulations, regardless of the control parameters $Ra_{d}$ and $H/d$. We
proceed to analyse this more quantitatively in the following section.

##  5 Mixing length and concentration profiles

The mixing length is defined as the vertical extent of the region where solute
concentration is non-uniform. Multiple ways to quantify the mixing length have
been proposed (Cook et0al., 2004), depending on whether a bulk-focused measure
is taken or it tracks the spread of the fastest growing fingers (e.g.,
$0.01\leq\overline{C}/C_{0}\leq 0.99$, see Gopalakrishnan et0al., 2017). The
latter method is used here to determine the mixing length obtained from the
experiments, where a threshold value corresponding to $0.96$ is applied to the
horizontally-averaged corrected light intensity as defined in0(13),
$\overline{I}_{2}$. The mixing length is then computed assuming that the flow
is symmetric with respect to the domain centerline, and a time correction is
also applied (see0§A.3).

Figure 8:  Dimensionless mixing length scaled with respect to $\ell$ for
experiments (colorbar) and simulations (legend). Simulations are shown for the
regular pattern (solid lines) and the randomly perturbed pattern (dotted
lines) as detailed in figure 4. In the experiments, mixing length is computed
assuming that the flow is symmetric with respect to the domain centerline, and
a time correction is also applied (see0§A.3). Data is only plotted up to the
time when the fingers reach the edge of the domain, so when $h=H$.
Asymptotically, experiments and simulations follow the scaling $h=2Ut$ (dashed
line).

The threshold-based definition discussed above would track the vertical
extremes of the finger growth over time. Alternatively, we can define a mixing
length based on the mean concentration profile, and this is the approach we
use for the simulations. For this, we can assume that the mean (i.e,
horizontally-averaged) profile takes a piecewise linear profile

$\frac{\overline{C}(z,t)}{C_{0}}=\begin{cases}0&z\leq-h/2\\\
1/2+z/h&|z|<h/2\\\ 1&z\geq h/2\end{cases}.$ (17)

Here, the overbar denotes a horizontal average, $h$ is the mixing length, and
the initial interface position is taken as $z=0$. This assumed profile
satisfies the integral relation

$h=\frac{6}{C_{0}^{2}}\int_{-\infty}^{\infty}\overline{C}(1-\overline{C})\,\mathrm{d}z,$
(18)

which we can use as a systematic definition of the mixing length $h$ in the
simulations. This integral definition provides a more useful statistical
measure of the mixing layer than the threshold-based definition, but is
impossible to compute with good accuracy from the experiments due to the
considerable thickness of the cell and the opacity of the dye.

The dimensionless mixing length scaled with respect to $\ell$ is shown for
experiments (colorbar) and simulations (legend) in figure08. Simulations are
presented for both the regular bead pattern (solid lines) and the perturbed
bead pattern (dotted lines) as shown in figure04. Asymptotically, both
experiments and simulations follow the scaling $h=2Ut$ (dashed line). The
buoyancy velocity $U$ represents the terminal velocity of a rising (falling)
parcel of light (heavy) fluid surrounded by heavy (light) fluid, and it is
achieved at the equilibrium between the driving force, represented by
buoyancy, and the drag due to viscous forces within the medium. This model is
a simplified representation of the flow, and any diffusion effect or
interaction with the flow structures is neglected. However, it provides a good
first order estimate for the growth of the mixing region. The observed
behaviour gives a precise indication of the flow dynamics: after a starting
phase in which the flow is influenced by the initial condition, the late-stage
dynamics is controlled by convection, with the mixing layer growing at the
characteristic buoyancy velocity $U$. This result is in agreement with the
observations in two-dimensional and three-dimensional Darcy simulations at
large Rayleigh-Darcy numbers (Borgnino et0al., 2021). At moderate Rayleigh-
Darcy numbers, a super-linear scaling has been observed (De0Paoli et0al.,
2019b), which has been interpreted has a finite size effect (Boffetta &
Musacchio, 2022; De0Paoli et0al., 2022a), i.e., the time/space available for
the fingers before touching the horizontal boundaries is insufficient to reach
this asymptotic regime.

Figure 9:  Vertical profiles of the horizontally-averaged concentration field
$\overline{C}(z,t)$ for simulation S12: $(a)$0plotted over the domain height;
$(b)$0plotted against the diffusive similarity variable; $(c)$0plotted against
height rescaled using the buoyancy velocity.

In the experiments, variations of solute concentration within the mixing layer
cannot be inferred with good accuracy. By contrast in the simulations, we have
all the information to quantify how the solute is distributed across the
mixing region. Specifically, in figure 9 we show the time evolution of the
horizontally averaged concentration profile $\overline{C}(z,t)$ for a given
simulation, in this case S12. This averaging is only performed over the fluid
fraction of the domain. In figure 9(a), we present the mean concentration over
the full height of the domain. The profile is constant in the regions above
and below the mixing layer, and transitions between 0 and $C_{0}$ in between,
with this mixing layer expanding over time. At very early times, before the
emergence of the Rayleigh-Taylor instability, we expect the mean profile to
develop diffusively. This is confirmed by figure 9(b), where the early-time
profiles of mean concentration are plotted against a rescaled spatial variable
$z/\sqrt{Dt}$. The profiles collapse, suggesting a self-similar development at
this stage, and agree well with the analytic solution for a diffusing
interface shown by the dashed black line. The profiles do not perfectly tend
to 0 and $C_{0}$ in this panel since the thin diffusive interface is contained
within the region seeded with initial noise. Once the Rayleigh-Taylor
instability develops and saturates, the dynamics are controlled by a balance
between the buoyancy driving and the drag provided by the porous medium. We
therefore expect the buoyancy velocity scale $U$ as defined in (4) to play an
important role in the spread of the solute. Indeed, by plotting the mean
concentration against the rescaled dimensionless coordinate $z/Ut$, we observe
further a self-similar behaviour, with $\overline{C}$ remaining close to a
linear profile within the mixing layer. Since the result of figure09(c) is
consistent with the assumption of (17), we are confident that the resulting
mixing length expression (18) is reliable for our simulations.

##  6 Flow structure and wavenumber

The flow structure is first discussed qualitatively with the aid of
experimental measurements, and then quantitatively using the numerical
results.

Figure 10:  Interface and flow structure evolution for two different
experiments. The instantaneous fluid-fluid interface is reported from early
(purple) to late (brown) times. For visualisation purposes, only a portion of
the entire domain width is shown. Experiment E2 is shown in (a), where the
inset corresponds to a squared region of side $H/20$. Panels0(c) and (d)
report the space-time intensity maps taken at $z/H=-0.05$ and $z/H=-0.10$,
respectively, as indicated with the black arrows in panel0(a), obtained from
the local values of raw light intensity ($I$) at that height $z$ normalised by
the maximum value within the picture ($I_{\text{max}}$). The time-dependent
interface evolution for experiment E12 is reported in panel0(d), where the
inset corresponds to squared region of side $H/10$. Here the width of the flow
structures is comparable with the interstitial space, i.e. with the diameter
of the beads. Space-time maps of normalised intensity taken at different
heights (panels0e,f) reveal more clearly that in this case the flow has a
behaviour that is very different from E2, with thin fingers penetrating
individually into the unmixed region.

In figure010 we consider the interface evolution for two different
experiments, namely E2 and E12. In figure010(a), the evolution of the fluid-
fluid interface is reported from early (blue) to late (brown) times. The
interface is computed from the normalised intensity fields $I_{1}$, as
discussed in §3.1. The flow around the beads is recorded at a high resolution,
and the position of the interface is reconstructed. An example of the
behaviour of the interface at the scale of the pores is reported in the inset
of figure010(a), where a squared region with side $H/20$ is magnified. From
this qualitative view of the evolution of the interface, we observe that the
structures of the convective flow are much larger than the beads diameter. To
more quantitatively evaluate the interface shape at different heights, we
consider the space-time maps in figures010(b) and (c), taken at $z/H=-0.05$
and $z/H=-0.10$, respectively, as indicated with the black arrows in
figure010(a). These maps are obtained from the local values of raw light
intensity ($I$) at that height $z$ normalised by the maximum value within the
picture ($I_{\text{max}}$). We observe that, at both heights, the
characteristic wavelength of the interface is larger than the beads size,
which is a feature typically observed in Darcy-type flows. In addition, we
also observe that the lower the measurement location, the later the interface
appears, and the larger the amplitude of the interfacial oscillations.
Similarly, the time-dependent interface evolution for experiment E12 is
reported in figure010(d). In this case, the flow corresponds to a
high-$\operatorname{\mathit{Pe}}$ flow, and the wavelength of the interfacial
structures are comparable in size with the interstitial space, i.e. with the
diameter of the beads. Space-time maps of normalised intensity taken at
different heights (figures010e,f) reveal more clearly that in this case the
flow has a behaviour that is very different from the previous one, with thin
fingers penetrating into the unmixed region. This analysis, although
qualitative, provides information about the flow structure in different
conditions, from a dissipation-controlled flow (E2) to a buoyancy-controlled
flow (E12).

Figure 11: (a) Space-time plot of the concentration profile $C(x,t)$ just
above the initial interface height $z=0$ from simulation S12. (b) Time series
of the mid-height finger wavelength as calculated using (19) for each
simulation with a regular bead pattern.

We now turn to the simulations to provide some quantitative results on the
wavelength of the finger structures. Since we can only capture the lower edge
of the interface in the experiments, it is not possible to investigate the
evolving dynamics of the fingers at the centre of the mixing layer. By
contrast, in the simulations we have the full details of the flow field, and
can make more quantitative statements about the time evolution of the finger
structures. We track the evolution of the finger width by considering the
concentration profile in horizontal cuts near the initial interface position
$z=0$. In the cases where we position the solid beads in a regular hexagonal
packing, we can take a horizontal cut through the domain that only contains
the liquid phase. This is important for analysis of the finger width, for
which we will take advantage of the horizontal periodicity and use Fourier
transforms. Specifically, we take the Fourier transform of the concentration
field in the $x$-direction at a fixed height $z_{m}$ to obtain $\hat{C}(k,t)$,
and then define the mean horizontal wavelength of the plumes as

$\lambda(t)=\frac{\int_{0}^{\infty}k^{-1}\left|\hat{C}(k,t)\right|^{2}\,\mathrm{d}k}{\int_{0}^{\infty}\left|\hat{C}(k,t)\right|^{2}\,\mathrm{d}k}.$
(19)

As shown in figure 11, the finger structures at the centre of the domain
exhibit a coarsening behaviour, with the mean wavelength increasing over time.
This coarsening is consistent across all the simulations at varying $Ra_{d}$
and $H/d$. After initial transient phase related to the onset of the Rayleigh-
Taylor instability, the wavelength exhibits a power-law scaling of $t^{1/2}$.
Such a scaling is often associated with a diffusive process, although in
figure011 the collapse is observed in terms of the advective scale $U$ rather
than the diffusivity $D$. This result of a scaling close to $t^{1/2}$ is in
fair agreement with the Darcy simulations of De0Paoli et0al. (2019a), where
the mean wavenumber $\kappa\sim 1/\lambda$ exhibits a $t^{-0.6}$ scaling
during the growth of the mixing layer. Boffetta et0al. (2020) also observe
$t^{1/2}$ scaling for the horizontal wavenumber in 3D Darcy simulations of
Rayleigh-Taylor instability, finding a collapse with $\lambda\sim\sqrt{Dt}$.
This contrasts to our pore-scale simulations, where the rate of coarsening is
independent of the molecular diffusivity. As the mixing layer grows in our
simulations, and the fingers coarsen to scales larger than a few multiples of
the pore scale, the Darcy assumption seems to become more relevant, as would
be expected.

##  7 Mean scalar dissipation

Numerical simulations allow accurate local measurements of concentration
gradients. They are used here to infer the local mixing state of the system
via the mean scalar dissipation,

(20) χ=D⟨|∇C|2⟩,

where $\langle\cdot\rangle$ stands for volume average over the fluid volume.
The evolution observed for $\chi$ is similar across all the considered
simulations, and it is illustrated in figure012 for simulation S6. We observe
that four main regimes can be identified here: (i) an initial diffusive
regime, followed by (ii) a linear instability, (iii) a convection-driven
regime, and (iv) a final shutdown phase. In this section, we will discuss in
detail these phases of the mixing process.

Figure 12:  Evolution of the mean scalar dissipation ($\chi$) over time ($t$)
for simulation S6. The time is made dimensionless with $(d/U)$. Three
different simulations are shown, corresponding to the parameters as reported
in table03.2 and different arrangements of the obstacles (regular and
perturbed patterns, corresponding to solid and dotted lines respectively). A
few concentration fields in different phases of the process are also shown
(regular pattern).

The mean scalar dissipation can be related to other global quantities of the
flow (Hidalgo et0al., 2012; De0Paoli, 2023), such as the dissolution rate (in
case of systems permeable to solute at the boundaries) or the volume-averaged
squared concentration, $\langle C^{2}\rangle$. Given our no flux boundary
conditions for concentration, the relation

$\partial_{t}\langle C^{2}\rangle=-2\chi,$ (21)

holds and will be used in the following to describe the evolution of the
dissipation rate.

We now consider the volume-averaged scalar dissipation rate over time for all
simulations, reported in figure013. A natural length scale to be used to
analyse the results during the diffusive regime is $\ell$, defined in (5), and
therefore time is scaled with $\phi\ell/U$. In this frame, all simulations
nicely collapse onto the same curve in the initial diffusive phase, and we
provide in the following an analytical description of the mixing process.

Figure 13:  Dimensionless mean scalar dissipation rate over dimensionless time
obtained from the numerical simulations (with regular and perturbed bead
patterns, corresponding to solid and dotted lines respectively). Black dashed
line indicates the diffusive evolution derived in equation0(24).

The fluid is initially motionless, with a step-like concentration field. As a
result, one can neglect any velocity contribution and assume that the
concentration field is uniform in horizontal direction. It follows that the
initial development, $t/(\phi\ell/U)<3\times 10^{3}$, is described by a purely
diffusive 1-D solution, where

$C=\frac{C_{0}}{2}\left[1+\mathrm{erf}\left(\frac{z}{2\sqrt{Dt}}\right)\right].$
(22)

From this expression, we can derive the scalar dissipation rate by

$\chi=D\langle|\nabla
C|^{2}\rangle=\frac{D}{H}\int_{-H/2}^{H/2}|\partial_{z}C|^{2}\mathrm{d}z\approx\frac{D}{H}\int_{-\infty}^{\infty}|\partial_{z}C|^{2}\mathrm{d}z=\sqrt{\frac{D}{8\pi
t}}\frac{C_{0}^{2}}{H}.$ (23)

This expression made dimensionless with $UC_{0}^{2}/(\phi H)$ can be expressed
as a function of $t/(\phi\ell/U)$ as:

$\frac{\phi\chi}{UC_{0}^{2}/H}=\frac{1}{\sqrt{8\pi}}\left(\frac{t}{\phi\ell/U}\right)^{-1/2}.$
(24)

The corresponding behaviour is shown in figure013 (black dashed line). It is
in agreement with the numerical findings since all simulations, regardless of
their $Ra_{d}$ and $H/d$, nicely collapse onto the analytical prediction. The
same solution obtained in equation0(23) can be derived using equation0(21)
with the diffusive solution0(22), which gives:

$\chi=\sqrt{\frac{D}{8\pi t}}\frac{C_{0}^{2}}{H}g(\xi),\quad\text{with}\quad
g(\xi)=\text{erf}(\xi)-\text{erf}\left(\frac{\xi}{\sqrt{2}}\right)\sqrt{2e^{-\xi^{2}}}$
(25)

and $\xi=H/(2\sqrt{Dt})$. Note that during the diffusive regime $\xi\gg 1$ and
$g(\xi)\approx 1$, consistent with (23).

The system leaves this diffusive behaviour as soon as convective instabilities
develop, i.e., when finger-like structures form. These structures stretch the
interface, providing a larger area for diffusion to act over and thus
promoting mixing (with a corresponding increase in $\chi$). The time at which
the instabilities become significant depends both on the magnitude of the
initial perturbation applied to the concentration field, and on the extension
of the region at which this perturbation is applied.

At later times, buoyancy-driven fingers formed at the initial fluid-fluid
interface grow towards the horizontal walls. This process is on one hand
controlled by convection, which drives the fingers elongation and the
corresponding increase of their interfacial extension, and on the other hand
by diffusion, which reduces the concentration gradient across the fingers
interface in time (Gopalakrishnan et0al., 2017). In addition, the presence of
solid obstacles makes the fingers more prone to split, possibly increasing
further the fluid mixing.

We provide a first estimate for the maximum dissipation rate in the convective
regime by considering that dissipation mainly takes place within the mixing
layer, where local non-zero concentration gradients exist. Outside the mixing
layer the fluid is nearly homogeneous in concentration. The volume-averaged
dissipation rate can therefore be approximated as

$\chi=D\langle|\nabla C|^{2}\rangle\approx D\frac{h}{H}\langle|\nabla
C|^{2}\rangle_{ML},$ (26)

where $\langle\cdot\rangle_{ML}$ denotes an average value across the mixing
layer. By assuming that the convective fingers stretch the interface around
the mixing layer, we can approximate the mean scalar gradient in the mixing
layer as the gradient of (22) at the interface, thus as

$|\nabla C|\approx\frac{C_{0}}{2\sqrt{\pi Dt}}.$ (27)

Combining these approximations with the result of §5 for the growth of the
mixing region $(h\approx 2Ut)$, we arrive at the estimate

$\chi\approx D\frac{2Ut}{H}\frac{C_{0}^{2}}{4\pi
Dt}=\frac{1}{2\pi}\frac{UC_{0}^{2}}{H}.$ (28)

This expression (black dotted line) proves to be an overestimate in figure014.
The overestimation is expected since, in the frame of the interface, the
approximation given by (27) is the _maximum_ gradient rather than the average
over a certain scale.

Similarly to what is observed during the diffusive regime, the dissipation
rate can also be estimated here from the bulk squared concentration. In the
following argument, we take the exact relation (21) and assume that
$d_{t}\langle C^{2}\rangle\approx d_{t}\langle\overline{C}^{2}\rangle$. This
implies that the leading order effect of dissipation is on mixing the mean
concentration $\overline{C}$, but does _not_ mean that dominant contribution
to the dissipation is _from_ the mean profile. We recall from figure09(c)
that, during the convective regime, the mean concentration field is linear
within the mixing region. Then we can approximate the horizontally-averaged
concentration as:

$\frac{\overline{C}(z,t)}{C_{0}}=\begin{cases}0&\text{for}-H/2\leq
z\leq-h/2\\\ 1/2+z/(2Ut)&\text{for }|z|<h/2\\\ 1&\text{for }h/2\leq z\leq
H/2\end{cases},$ (29)

which is nothing but (17) with $h=2Ut$. Using equation0(21) with
approximation0(29), we estimate the mean scalar dissipation as

$\chi\approx-\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{-H/2}^{-h(t)/2}\overline{C}^{2}\mathrm{d}z+\int_{-h(t)/2}^{h(t)/2}\overline{C}^{2}\mathrm{d}z+\int_{h(t)/2}^{H/2}\overline{C}^{2}\mathrm{d}z\right)=\frac{1}{6}\frac{UC_{0}^{2}}{H}.$
(30)

This result, shown as a red dotted line in figure014, is analogous to the
expression derived in equation0(28). The obtained numerical coefficient is
also similar. As mentioned above, this estimate does not take into account
local fluctuations of concentration, likely responsible of the decrease
observed for the scalar dissipation, and to that aim higher order statistics
should be considered.

Figure 14:  Dimensionless mean scalar dissipation rate over dimensionless time
obtained from the numerical simulations (with regular and perturbed bead
patterns, corresponding to solid and dotted lines respectively, colors as in
Fig.013). The black dotted line, corresponding to a dimensionless value of
$1/2\pi$ from (28), refers to the maximum mean dissipation rate in the
convective regime. The alternative estimate of $1/6$ derived in (30) is also
shown as a red dotted line. The thicker black dashed line denotes the time-
dependent decreasing mean dissipation, as predicted by (34), in the same
regime. The fitting parameter $\beta=0.75$ is chosen such that the estimate
tends to 0.1 as $t\rightarrow\infty$.

When plotted on a linear scale, as in figure014, it becomes clear that the
dissipation rate decreases with time. The mechanism that controls the
dissipation rate in this phase is possibly due to several interplaying
phenomena. While the mixing region grows at a constant speed ($h\approx 2Ut$),
the concentration gradients across the fingers’ interface reduce, owing to
diffusion. In addition, the number of fingers varies in time, as well as the
extension of the interfacial region across which the fluids mix. We propose
here a model that takes into account these features of the flow and explains
the decreasing behaviour of $\chi$ in the convective phase.

A schematic model for the fluid-fluid interface is proposed figure015(a). For
any instant $t$ in the convective regime, fingers are considered to have all
the same length, which matches the values of the mixing length $h(t)$. In
addition, we also consider the fingers to have all the same width, $\lambda$.
We approximate the fingers shape to be straight, and as a result the length of
the interface $L_{i}$ to be equal to the segmented blue line in figure015(a),
which is given by:

$L_{i}(t)=L+2\frac{L}{\lambda(t)}h(t).$ (31)

We introduce in figure015(b) a coordinate system defined such that $s$ is
tangential to the interface and $n$ perpendicular to it. As a result, assuming
an interface with uniform thickness $\delta(t)$, the mean scalar dissipation
can be computed by integrating the $|\partial_{n}C|^{2}$ across the thickness
$\delta$ and along the interface length $L_{i}$ to obtain

$\chi=D\langle|\nabla
C|^{2}\rangle=\frac{DL_{i}}{HL}\int_{-\delta/2}^{+\delta/2}|\partial_{n}C|^{2}\mathop{}\\!\mathrm{d}n\approx\frac{DL_{i}}{HL}\frac{C_{0}^{2}\delta}{4\pi
Dt}.$ (32)

Here, we assumed that the concentration gradient across the interface evolves
diffusively using (27), and took $\delta/(2\sqrt{Dt})\ll 1$. Motivated by the
observations of coarsening fingers in figure 11(b), we assume a diffusive
growth for $\delta$ with an effective diffusivity of $dU$, namely:

$\delta=\beta\sqrt{dUt},$ (33)

with $\beta$ being a fitting parameter. Note that the effective diffusivity
$dU$ determined here coincides with the effective longitudinal dispersion
usually considered for bead packs (Tsinober et0al., 2023; Liang et0al., 2018).
We now substitute the expressions obtained in (31) and (33) into the
definition (32), together with the evolution of the mixing length $h=2Ut$ and
the wavelength measured from the simulations $\lambda/d=\alpha\sqrt{t/(d/U)}$,
as in figure011. Finally, we obtain an expression for $\chi(t)$ in the
convective regime:

$\chi=\frac{C_{0}^{2}U}{H}\frac{\beta}{\alpha\pi}\left[1+\frac{\alpha}{4}\left(\frac{t}{d/U}\right)^{-1/2}\right],$
(34)

in which $\alpha=2.39$ is obtained from the numerical results of figure 11 and
$\beta=0.75$ is obtained as a fitting parameter for the data of $\chi(t)$.

Figure 15:  Model for the evolution of the interface in the convective phase.
(a)0The entire domain (dimensions $L\times H$) is sketched for a time $t$
within the convective regime, when fingers have already developed. The average
width of the fingers is $\lambda$. At this time, the extension of the mixing
region is $h(t)$, and the interface between the fluids (red line) can be
approximated by a segmented line (blue). (b)0Detail of the interfacial region.
We assume the interface (black solid line) has a finite and uniform thickness
$\delta(t)$, and the mixing occurs within this region. A coordinate system is
defined such that $s$ is tangential to the interface and $n$ perpendicular to
it.

Comparison of the analytical prediction (34) (thick black dashed line) against
the numerical results of mean scalar dissipation is shown in figure014. This
solution captures more accurately the evolution of the mixing process, and it
is also quantitatively in agreement with the maximum dissipation (28) (thin
black dashed line) at early times. This result suggests also that mixing
process during the convective regime, unlike the diffusive regime, is
controlled by the characteristic pore-scale of the domain, $d$. The solution
is indeed independent of the domain size, $H/d$, and the driving force,
$Ra_{d}$.

We observe that the behaviour of $\chi$ shown in figure014 is remarkably
dissimilar compared to the corresponding Darcy simulations (De0Paoli et0al.,
2019a), where it is observed to increase with time up to impingement of the
fingers on the horizontal boundaries. Hidalgo et0al. (2015) observed that in
presence of fluids characterized by non-monotonic density-concentration
profiles, which leads to a considerably different flow topology, the scalar
dissipation remains constant during the convective regime. We speculate that
the pore-induced dispersion effects (Dentz et0al., 2018) are likely
responsible for these differences.

Figure 16:  (a) Dimensionless mean scalar dissipation rate against
dimensionless time obtained from numerical simulations (with regular and
perturbed bead patterns, corresponding to solid and dotted lines respectively,
colors as in Fig.013). The horizontal dashed line, corresponding to a
dimensionless value of $1/2\pi$, equation0(28), refers to the maximum mean
dissipation rate in the convective regime. Vertical dashed lines approximately
correspond to the time required for the fingers to reach the horizontal walls
of the domain, $t/(H/U)=1/2$, and for the core of the domain to be influenced
by the presence of the walls, $t/(H/U)=1$. (b) Vertical profiles of the root-
mean-squared concentration from simulation S12, highlighting the uniform
region of concentration variance within the region $|z|<Ut/2$. Instantaneous
profiles are plotted for $t/(U/d)\geq 5$, each with an opacity of $0.1$.

During the convective regime, the fingers grow under the action of buoyancy
and eventually touch the horizontal boundaries of the domain. This event
occurs approximately at ${t=H/(2U)}$, given the growth rate of the mixing
region, and marks a first observable reduction of the scalar dissipation. It
is apparent that in this regime, referred here to as a shutdown regime, the
relevant flow length scale is $H$. Therefore we report in figure016(a) the
evolution of $\chi$ as a function of $t/(H/U)$, and we observe that all curves
nicely collapse in the late stage of the flow evolution, i.e. for $0.25\leq
t/(H/U)\leq 1$.

After the fingers impinge on the horizontal walls, the concentration field in
the near-wall regions begins to be progressively more homogeneous. This
reflects the reduction of the local concentration gradients, and therefore of
the mean scalar dissipation. Mixing is still ongoing, however, in the central
portion of the domain, where the information that the walls are present has
not reached yet. This can be seen in figure 16(b), where we plot the root-
mean-squared concentration profiles for a typical simulation. With height
rescaled by $Ut$, we observe a ‘core’ region of uniform scalar variance
persists between $-Ut/2\leq z\leq Ut/2$. The shape of this profile persists
even after the fingers reach the domain boundaries. Approximately at time
$t=H/U$, the core of the flow is affected as well, and the entire domain
becomes more homogeneous in solute concentration, i.e., the local
concentration gradients are small and the mean scalar dissipation drops
further. The overall dynamics is controlled in this case by the domain height,
however geometry and buoyancy still play a role in determining the reduction
rate of the scalar dissipation.

##  8 Summary, conclusions and outlook

We have studied the process of convective dissolution in porous media. We
considered a Rayleigh-Taylor instability, consisting of two miscible fluid
layers of different density placed in an unstable configuration, with the
heavy fluid on top of the lighter one. The flow is unstable due to the
presence of a solute, which induces the density differences driving the mixing
process. The porous medium consists of a confined, homogeneous and isotropic
bead pack, with the solid obstacles being impermeable to fluid and solute. We
investigated the flow at the scale of the pores using experimental
measurements, numerical simulations and physical models. Simulations employ
finite differences coupled to the immersed boundary method, while experiments
are performed with transparent beads and fully miscible fluids. Experiments
and simulations have been specifically designed to mimic the same flow
conditions, namely linear dependency of fluid density with solute
concentration, high Schmidt numbers ($\operatorname{\mathit{Sc}}=O(10^{2})$)
and the same values of porosity.

The evolution of the flow is quantified by the mixing length $h$, which
represents the vertical extension of the mixing region. We demonstrate via
experiments and simulations that the system is characterized by a linear
scaling $h=2Ut$, and the growth of the mixing region is controlled by the
buoyancy velocity $U$. This velocity is achieved at the equilibrium between
driving forces (density contrast between the fluids, $\Delta\rho$) and viscous
dissipation (drag through the medium). The solute evolution presents a self-
similar behaviour during the initial diffusive phase and during the following
convective regime. We demonstrate this self-similarity of the flow by
inspection of the horizontally-averaged concentration profiles. The flow
structures are analysed using the mean wavelength at the centreline
($\lambda$), and also in this case a scaling holds ($\lambda\sim\sqrt{Ut}$).
Finally, we analyse the mixing dynamics of the system, which is quantified via
the mean scalar dissipation rate, $\chi$, and three flow regimes are observed.
The flow is initially controlled by diffusion $(\chi\sim\sqrt{D/t})$, which
produces solute mixing across the initial horizontal interface. When the
interfacial diffusive layer is sufficiently thick and it eventually becomes
unstable, finger-like structures form and drive the flow into a convection-
dominated phase. In this case, fluid mixing is controlled by diffusion
(predominantly at the sides of the fingers) and by buoyancy-driven convection
(which makes the finger grow vertically as $h=2Ut$). After an initial growth
at the end of which the dissipation rate attains a maximum value (which we
predict), a reduction of the mean dissipation rate is observed as a result of
the competing effect of merging of the fingers (negative contribution) that
dominates over the growth of the mixing region (positive contribution). To
describe this behaviour, we propose a model of mixing that relies on the
diffusion process across the interface of the fingers. Finally, when the
fingers grow sufficiently to touch the horizontal boundaries of the domain,
the scalar dissipation reduces dramatically (shutdown regime), due to the
absence of fresh fluid contributing to mixing. We further elucidated the
physics of the observed phenomena with the aid of simple physical models, and
we demonstrated that each regime is controlled by a different length scale,
namely the scale of diffusion, the characteristic length of the pores, or the
domain height.

In this study we have focused on modelling the mixing dynamics and the flow
structure of porous media convection. We used simulations performed in two-
dimensional arrays of circular objects and experiments consisting of three-
dimensional packs of spherical beads. In a porous Rayleigh-Taylor system, at
the Darcy scale, the mixing region is observed to grow faster in two-
dimensional domains than in three-dimensional ones (Borgnino et0al., 2021).
Unlike in the turbulent case, the transition between these regimes occurs
sharply when the thickness of the domain exceeds the wavelength of the most
unstable mode (Boffetta & Musacchio, 2022). The nature of this transition in
pore-resolved flows has not been explored yet, and in the future it will be of
interest to extend the present study to simulations in three dimensions, as
well as two-dimensional experiments (De et0al., 2022), to allow for a one-to-
one comparison of these findings and to investigate the effect of the
dimensionality of the flow on the evolution of a pore-resolved system.

Present findings are relative to domains having a dimension up to few hundred
pores. In contrast, geophysical applications involve formations that may be
orders of magnitude larger (Huppert & Neufeld, 2014), and modelling of the
entire dissolution process may be required (Wang et0al., 2022; Szulczewski
et0al., 2013). For instance, in the case of carbon sequestration it is desired
to determine the time required to dissolved a prescribed fraction of the
injected fluid volume, or to estimate the horizontal spread of the current of
CO2 (MacMinn et0al., 2012; De0Paoli, 2021). Performing pore-scale studies of
such flows is beyond the present capabilities and Darcy simulations including
the effect of dispersion are required (Dentz et0al., 2023). In this context,
one remarkable finding of the present study is the behaviour of the mean
scalar dissipation during the convection-dominated regime, which is observed
to decrease in time in pore-resolved flows (see figure014) while it grows in
Darcy flows (De0Paoli et0al., 2019a). A different behaviour is observed by
Hidalgo et0al. (2015) in the presence of fluids characterized by non-monotonic
density-concentration profiles. In their Darcy simulations, this property of
the fluids leads to a remarkably different flow topology with a nearly flat
interface between the two fluid layers, and the scalar dissipation remains
constant during the convective regime. We speculate that the pore-induced
dispersion effects (Dentz et0al., 2018), not captured by Darcy simulations,
are likely responsible for these differences and should be accounted by
suitable dispersion models developed from pore-resolved data.

The results presented are relevant to homogeneous and isotropic porous media
at the porosity $\phi=0.37$. Geological formations are generally characterized
by a more complex structure, typically heterogeneous at multiple scales
(Hidalgo et0al., 2022), anisotropic (Ennis-King et0al., 2005) and with
porosity values that can be as low as 0.05 (Bickle et0al., 2017). Exploring
the effect of such irregularities of the medium on the flow structure would
represent an important extension of this study to be considered for future
works.

[Funding] The work of CJH was funded by the Max Planck Center for Complex
Fluid Dynamics. We acknowledge PRACE for awarding us access to MareNostrum4 at
the Barcelona Supercomputing Center (BSC), Spain and IRENE at Très Grand
Centre de Calcul (TGCC) du CEA, France (project 2021250115). This work was
also carried out on the Dutch national e-infrastructure with the support of
SURF Cooperative. [Declaration of interests]The authors report no conflict of
interest. [Author ORCID] Marco De Paoli,
https://orcid.org/0000-0002-4709-4185; Christopher J. Howland
https://orcid.org/0000-0003-3686-9253; Roberto Verzicco
https://orcid.org/0000-0002-2690-9998; Detlef Lohse
https://orcid.org/0000-0003-4138-2255.

##  Appendix A Additional experimental details

###  A.1 Fluid density

Following the empirical correlations of Novotný & Söhnel (1988), the density
of an aqueous solution of KMnO4 can be determined as a function of the solute
concentration $C$ (expressed in mol/dm3) and the temperature of the solution
$\vartheta$ (expressed in ${}^{\circ}\text{C}$) as:

$\rho(C,\vartheta)=\rho_{w}(\vartheta)+A_{1}(\vartheta)C+A_{2}(\vartheta)C^{3/2}\text{
,}$ (35)

where the water density $\rho_{w}(\vartheta)$ and the temperature-dependent
coefficients $A_{1}(\vartheta)$ and $A_{2}(\vartheta)$, are given by:

$\displaystyle\rho_{w}(\vartheta)$ $\displaystyle=$ $\displaystyle+9.997\times
10^{2}\frac{\text{kg}}{\text{m}^{3}}+2.044\times
10^{-1}\frac{\text{kg}}{{}^{\circ}\text{C}\cdot\text{m}^{3}}\cdot\vartheta+$
(36) $\displaystyle-6.174\times
10^{-2}\frac{\text{kg}}{(^{\circ}\text{C})^{3/2}\cdot\text{m}^{3}}\cdot\vartheta^{3/2}\text{
,}$ $\displaystyle A_{1}(\vartheta)$ $\displaystyle=$
$\displaystyle+1.223\times 10^{-1}\frac{\text{kg}}{\text{mol}}-1.029\times
10^{-4}\frac{\text{kg}}{{}^{\circ}\text{C}\text{ mol}}\cdot\vartheta+$ (37)
$\displaystyle+8.093\times
10^{-6}\frac{\text{kg}}{(^{\circ}\text{C})^{2}\text{
mol}}\cdot\vartheta^{2}\text{ ,}$ $\displaystyle A_{2}(\vartheta)$
$\displaystyle=$ $\displaystyle-1.485\times
10^{1}\frac{\text{kg}}{\text{m}^{3}}\left(\frac{\text{dm}^{3}}{\text{mol}}\right)^{3/2}+9.079\times
10^{-1}\frac{\text{kg}}{{}^{\circ}\text{C
}\text{m}^{3}}\left(\frac{\text{dm}^{3}}{\text{mol}}\right)^{3/2}\cdot\vartheta+$
(38) $\displaystyle-7.566\times 10^{-3}\frac{\text{kg}}{(^{\circ}\text{C
}^{2})\text{m}^{3}}\left(\frac{\text{dm}^{3}}{\text{mol}}\right)^{3/2}\cdot\vartheta^{2}\text{
.}$

The relative dependence of density of the solution, $\rho$, mass fraction,
$\omega$, and solute concentration, $C$, is shown in figure017 (blue solid
lines). We report here $C(\omega)$0(panel0a),
$\rho(\omega)-\rho_{w}$0(panel0b) and $\rho(C)-\rho_{w}$0(panel0c), being
$\rho_{w}$ the water density defined as in0(36). The density of the mixture
(figure017c), $\rho(C)$, is well approximated by a linear function of the
solute concentration (10) (dashed grey line).

Figure 17:  Given the definition0(9) and the empirical correlations0(35)-(38),
the relative dependence of density of the solution, $\rho$, mass fraction,
$\omega$, and solute concentration, $C$, is obtained (blue solid lines). We
report $C(\omega)$0(panel0a), $\rho(\omega)-\rho_{w}$0(panel0b) and
$\rho(C)-\rho_{w}$0(panel0c), being $\rho_{w}$ the water density defined as
in0(36). The profiles shown here correspond to $\vartheta=25^{\circ}\text{C}$.
The fluid properties are well approximated by linear functions (grey dashed
lines).

Figure 18: Beads characterization for all diameters considered. Left column
(i): distribution of the diameter of the beads normalized with respect to the
nominal diameter ($d_{n}$). Central column (ii): the eccentricity of the beads
($\varepsilon$) when approximated to ellipses. Right column: example of beads
employed. A reference length (red bar) corresponding to 50mm is also
indicated.

###  A.2 Characterization of the porous medium

The employed glass beads consist of clear, polished glass spheres manufactured
by various producers (Hecht Karl, Witeg). The size distribution of the beads
has been determined optically. The beads are placed in a transparent container
on top of a homogeneous light source (Phlox LEDW-BL $300\times 300$0mm2). The
images are recorded with a high-resolution camera (Nikon D 850 with lenses
Sigma DG Macro 1050mm). Details are reported in figure018. For all bead size
considered, the beads are measured by finding the best-fitting ellipse, and
the mean diameter is obtained as $d=\sqrt{ab}$, being $a,b$ the major and
minor axis of the ellipse, respectively. Finally, the shape is also evaluated
with the eccentricity, defined as $\varepsilon=a/b$. A summary of the beads
size and shape is reported in figure018. We observe from the statistics
(figures018c-i,c-ii) and from the pictures (figures018c-iii), that the beads
having nominal diameter $d_{n}=3$0mm present a wider distribution of diameters
and shapes compared to the other diameters (figures018a,b,d) .

###  A.3 Correction for the experimental measurement of the mixing length

The first instant considered ($t=0$) corresponds to first image analyzed after
the experiments has started, i.e. after (i)0the gate has been removed,
(ii)0the cell is turned upside-down, and (iii)0the cell is placed in front of
the camera. This operation produces a delay in the acquisition of the images,
and the original time at which the flow starts cannot be accurately
determined. In addition, the influence of the initial perturbation is
apparent. The initial interface may not be flat, and plumes may form due to
the gate opening or due to operational procedure, producing a perturbation of
the concentration field. Since we observe that after an initial transition the
evolution of the mixing length is linear, we correct the time used to report
the mixing length as follows. We first fit the mixing length using a linear
function, and then we determine the time $t_{0}$ at which the value of this
function is zero. The new time considered for the experiments is then obtained
by adding the quantity $t_{0}$ to the original dimensionless time
$t/(\phi\ell/U)$. A graphical representation of this correction process is
shown in figure019.

Figure 19:  Graphical representation of the correction procedure applied to
the experimental measurements. Experiments (solid line) are fitted by a linear
function (dashed line). The new time considered for the experiments is
obtained by adding the quantity $t_{0}$ to the original dimensionless time
$t/(\phi\ell/U)$.

## References

  * Ataei-Dadavi et0al. (2019) Ataei-Dadavi, I, Rounaghi, N, Chakkingal, M, Kenjeres, S, Kleijn, CR & Tummers, MJ 2019 An experimental study of flow and heat transfer in a differentially side heated cavity filled with coarse porous media. Int. J. Heat Mass Transf. 143, 118591\.
  * Backhaus et0al. (2011) Backhaus, S, Turitsyn, K & Ecke, RE 2011 Convective instability and mass transport of diffusion layers in a Hele-Shaw geometry. Phys. Rev. Lett. 1060(10), 104501\.
  * Bejan (2013) Bejan, A 2013 Convection heat transfer. John Wiley & Sons.
  * Bickle et0al. (2017) Bickle, M, Kampman, N, Chapman, H, Ballentine, C, Dubacq, B, Galy, A, Sirikitputtisak, T, Warr, O, Wigley, M & Zhou, Z 2017 Rapid reactions between CO2, brine and silicate minerals during geological carbon storage: Modelling based on a field CO2 injection experiment. Chem. Geol. 468, 17–31.
  * Boffetta et0al. (2020) Boffetta, G., Borgnino, M. & Musacchio, S. 2020 Scaling of Rayleigh-Taylor mixing in porous media. Phys. Rev. Fluids 50(6), 62501.
  * Boffetta & Mazzino (2017) Boffetta, G & Mazzino, A 2017 Incompressible Rayleigh-Taylor turbulence. Ann. Rev. Fluid Mech. 49, 119–143.
  * Boffetta & Musacchio (2022) Boffetta, G. & Musacchio, S. 2022 Dimensional effects in Rayleigh–Taylor mixing. Philos. Trans. Royal Soc. 3800(2219), 20210084.
  * Borgnino et0al. (2021) Borgnino, M, Boffetta, G & Musacchio, S 2021 Dimensional transition in Darcy-Rayleigh-Taylor mixing. Phys. Rev. Fluids 60(7), 074501.
  * Brouzet et0al. (2022) Brouzet, C., Méheust, Y. & Meunier, P. 2022 Co 2 convective dissolution in a three-dimensional granular porous medium: An experimental study. Phys. Rev. Fluids 70(3), 033802\.
  * Cardin & Olson (1994) Cardin, P & Olson, P 1994 Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core. Phys. Earth Planet. Inter. 820(3-4), 235–259.
  * Chakkingal et0al. (2019) Chakkingal, M, Kenjereš, S, Ataei-Dadavi, I, Tummers, MJ & Kleijn, CR 2019 Numerical analysis of natural convection with conjugate heat transfer in coarse-grained porous media. Int. J. Heat Fluid Fl. 77, 48–60.
  * Cook et0al. (2004) Cook, AW, Cabot, W & Miller, PL 2004 The mixing transition in Rayleigh–Taylor instability. J. Fluid Mech. 511, 333–362.
  * De et0al. (2022) De, N, Singh, N, Fulcrand, R, Méheust, Y, Meunier, P & Nadal, F 2022 Two-dimensional micromodels for studying the convective dissolution of carbon dioxide in 2D water-saturated porous media. Lab Chip 220(23), 4645–4655.
  * De0Paoli (2021) De0Paoli, M. 2021 Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide. Phys. Fluids 330(1), 016602.
  * De0Paoli (2023) De0Paoli, M. 2023 Convective mixing in porous media: A review of Darcy, pore-scale and Hele-Shaw studies, arXiv: 2310.01999.
  * De0Paoli et0al. (2020) De0Paoli, M, Alipour, M & Soldati, A 2020 How non-darcy effects influence scaling laws in Hele-Shaw convection experiments. J. Fluid Mech. 892, A41.
  * De0Paoli et0al. (2019a) De0Paoli, M, Giurgiu, V, Zonta, F & Soldati, A 2019a Universal behavior of scalar dissipation rate in confined porous media. Phys. Rev. Fluids 40(10), 101501\.
  * De0Paoli et0al. (2022a) De0Paoli, M., Perissutti, D., Marchioli, C. & Soldati, A. 2022a Experimental assessment of mixing layer scaling laws in Rayleigh-Taylor instability. Phys. Rev. Fluids 70(9), 093503.
  * De0Paoli et0al. (2022b) De0Paoli, M., Pirozzoli, S., Zonta, F. & Soldati, A. 2022b Strong Rayleigh–Darcy convection regime in three-dimensional porous media. J. Fluid Mech. 943, A51.
  * De0Paoli et0al. (2016) De0Paoli, M, Zonta, F & Soldati, A 2016 Influence of anisotropic permeability on convection in porous media: implications for geological CO2 sequestration. Phys. Fluids (1994-present) 280(5), 056601.
  * De0Paoli et0al. (2017) De0Paoli, M, Zonta, F & Soldati, A 2017 Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown. Phys Fluids (1994-present) 290(2), 026601\.
  * De0Paoli et0al. (2019b) De0Paoli, M, Zonta, F & Soldati, A 2019b Rayleigh-Taylor convective dissolution in confined porous media. Phys. Rev. Fluids 4, 023502.
  * De0Wit (2020) De0Wit, A 2020 Chemo-hydrodynamic patterns and instabilities. Ann. Rev. Fluid Mech. 52, 531–555.
  * Delgado (2007) Delgado, JMPQ 2007 Longitudinal and transverse dispersion in porous media. Chemical Engineering Research and Design 850(9), 1245–1252.
  * Dentz et0al. (2023) Dentz, M, Hidalgo, JJ & Lester, D 2023 Mixing in porous media: concepts and approaches across scales. Transp. Porous Media 1460(1-2), 5–53.
  * Dentz et0al. (2018) Dentz, M, Icardi, M & Hidalgo, JJ 2018 Mechanisms of dispersion in a porous medium. J. Fluid Mech. 841, 851–882.
  * Dhar et0al. (2022) Dhar, J., Meunier, P, Nadal, F & Méheust, Yves 2022 Convective dissolution of carbon dioxide in two-and three-dimensional porous media: The impact of hydrodynamic dispersion. Phys.Fluids 340(6), 064114.
  * Dullien (1991) Dullien, FAL 1991 Porous media: fluid transport and pore structure. Academic press.
  * Eckel et0al. (2022) Eckel, AM, Liyanage, R, Kurotori, T & Pini, R 2022 Spatial Moment Analysis of Convective Mixing in Three-Dimensional Porous Media Using X-ray CT Images. Ind. Eng. Chem. Res. .
  * Emami-Meybodi et0al. (2015) Emami-Meybodi, H., Hassanzadeh, H., Green, C.0P. & Ennis-King, J. 2015 Convective dissolution of CO2 in saline aquifers: Progress in modeling and experiments. Int. J. Greenh. Gas Con. 40, 238–266.
  * Ennis-King et0al. (2005) Ennis-King, J, Preston, I & Paterson, L 2005 Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys. Fluids (1994-present) 170(8), 084107\.
  * Fadlun et0al. (2000) Fadlun, E.0A., Verzicco, R., Orlandi, P. & Mohd-Yusof, J. 2000 Combined Immersed-Boundary Finite-Difference Methods for Three-Dimensional Complex Flow Simulations. Journal of Computational Physics 1610(1), 35–60.
  * Feltham et0al. (2006) Feltham, D.0L., Untersteiner, N., Wettlaufer, J.0S. & Worster, M.0G. 2006 Sea ice is a mushy layer. Geophys. Res. Lett. 330(14).
  * Fu et0al. (2013) Fu, X, Cueto-Felgueroso, L & Juanes, R 2013 Pattern formation and coarsening dynamics in three-dimensional convective mixing in porous media. Philos. T. R. Soc. A 3710(2004), 20120355\.
  * Gasow et0al. (2021) Gasow, S, Kuznetsov, AV, Avila, M & Jin, Y 2021 A macroscopic two-length-scale model for natural convection in porous media driven by a species-concentration gradient. J. Fluid Mech. 926, A8.
  * Gasow et0al. (2022) Gasow, S, Kuznetsov, AV & Jin, Y 2022 Prediction of pore-scale-property dependent natural convection in porous media at high rayleigh numbers. Int. J. Therm. Sci. 179, 107635\.
  * Gasow et0al. (2020) Gasow, S, Lin, Z, Zhang, HC, Kuznetsov, A0V, Avila, M & Jin, Y. 2020 Effects of pore scale on the macroscopic properties of natural convection in porous media. J. Fluid Mech. .
  * Gopalakrishnan et0al. (2017) Gopalakrishnan, S.0S., Carballido-Landeira, J., De0Wit, A. & Knaepen, B. 2017 Relative role of convective and diffusive mixing in the miscible Rayleigh-Taylor instability in porous media. Phys. Rev. Fluids 2, 012501.
  * Grossmann & Lohse (2000) Grossmann, S & Lohse, D 2000 Scaling in thermal convection: a unifying theory. J. Fluid Mech. 407, 27–56.
  * Grossmann & Lohse (2001) Grossmann, S & Lohse, D 2001 Thermal convection for large prandtl numbers. Phys. Rev. Lett. 860(15), 3316.
  * Happel & Brenner (2012) Happel, J. & Brenner, H. 2012 Low Reynolds number hydrodynamics: with special applications to particulate media, , vol.01\. Springer Science & Business Media.
  * Hartmann et0al. (2001) Hartmann, DL, Moy, LA & Fu, Q 2001 Tropical convection and the energy balance at the top of the atmosphere. J. Clim. 140(24), 4495–4511.
  * Haughey & Beveridge (1969) Haughey, DP & Beveridge, GSG 1969 Structural properties of packed beds - a review. Can. J. Chem. Eng. 470(2), 130–140.
  * Hewitt (2020) Hewitt, DR 2020 Vigorous convection in porous media. Proc. Math. Phys. Eng. Sci. 4760(2239), 20200111.
  * Hewitt et0al. (2012) Hewitt, DR, Neufeld, JA & Lister, JR 2012 Ultimate regime of high Rayleigh number convection in a porous medium. Phys. Rev. Lett. 1080(22), 224503.
  * Hewitt et0al. (2014) Hewitt, DR, Neufeld, JA & Lister, JR 2014 High Rayleigh number convection in a three-dimensional porous medium. J. Fluid Mech. 748, 879–895.
  * Hidalgo et0al. (2012) Hidalgo, JJ, Fe, J, Cueto-Felgueroso, L & Juanes, R 2012 Scaling of convective mixing in porous media. Phys. Rev. Lett. 1090(26), 264503.
  * Hidalgo et0al. (2022) Hidalgo, JJ, Neuweiler, I & Dentz, M 2022 Advective trapping in the flow through composite heterogeneous porous media. Transp. Porous Med. 1430(3), 599–618.
  * Hidalgo et0al. (2015) Hidalgo, J.0J., Dentz, M., Cabeza, Y. & Carrera, J. 2015 Dissolution patterns and mixing dynamics in unstable reactive flow. Geophys. Res. Lett. 420(15), 6357–6364.
  * Horton & Rogers0Jr (1945) Horton, CW & Rogers0Jr, FT 1945 Convection currents in a porous medium. Journal of Applied Physics 160(6), 367–370.
  * Huppert & Neufeld (2014) Huppert, HE & Neufeld, JA 2014 The fluid mechanics of carbon dioxide sequestration. Ann. Rev. Fluid Mech. 46, 255–272.
  * Korba & Li (2022) Korba, D & Li, L 2022 Effects of pore scale and conjugate heat transfer on thermal convection in porous media. J. Fluid Mech. 944, A28.
  * Krepper et0al. (2002) Krepper, E, Hicken, E-F & Jaegers, H 2002 Investigation of natural convection in large pools. Int. J. Heat Fluid Flow 230(3), 359–365.
  * Krevor et0al. (2023) Krevor, S, de0Coninck, H, Gasda, S0E, Ghaleigh, NS, de0Gooyert, V, Hajibeygi, H, Juanes, R, Neufeld, J, Roberts, J0J & Swennenhuis, F 2023 Subsurface carbon dioxide and hydrogen storage for a sustainable energy future. Nat. Rev. Earth Environ. pp. 1–17.
  * Lapwood (1948) Lapwood, ER 1948 Convection of a fluid in a porous medium. Math. Proc. Camb. Philos. Soc. 44, 508–521.
  * Lasser et0al. (2021) Lasser, J., Ernst, M. & Goehring, L. 2021 Stability and dynamics of convection in dry salt lakes. J. Fluid Mech. 917, A14.
  * Lasser et0al. (2023) Lasser, J., Nield, J.0M, Ernst, M., Karius, V., Wiggs, GFS, Threadgold, MR, Beaume, C & Goehring, L 2023 Salt polygons and porous media convection. Phys. Rev. X 130(1), 011025.
  * LeBlanc (1984) LeBlanc, DR 1984 Sewage plume in a sand and gravel aquifer, Cape Cod, Massachusetts, , vol. 2218. US Geological Survey.
  * Liang et0al. (2018) Liang, Y., Wen, B., Hesse, M.0A. & DiCarlo, D. 2018 Effect of dispersion on solutal convection in porous media. Geophys. Res. Lett. 450(18), 9690–9698.
  * Liu et0al. (2020) Liu, S, Jiang, L, Chong, KL, Zhu, X, Wan, Z0Hua, Verzicco, R, Stevens, R.J.A.M., Lohse, D & Sun, C 2020 From Rayleigh-Bénard convection to porous-media convection: How porosity affects heat transfer and flow structure. J. Fluid Mech. 895.
  * MacMinn et0al. (2012) MacMinn, CW, Neufeld, JA, Hesse, MA & Huppert, HE 2012 Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers. Water Resour. Res. 480(11).
  * Marshall & Schott (1999) Marshall, J & Schott, F 1999 Open-ocean convection: Observations, theory, and models. Rev. Geophys. 370(1), 1–64.
  * Menand & Woods (2005) Menand, T & Woods, AW 2005 Dispersion, scale, and time dependence of mixing zones under gravitationally stable and unstable displacements in porous media. Water Resourc. Res. 410(5).
  * Middleton et0al. (2016) Middleton, C.0A., Thomas, C., de0Wit, A. & Tison, J.-L. 2016 Visualizing brine channel development and convective processes during artificial sea-ice growth using Schlieren optical methods. J. Glaciol. 620(231), 1–17.
  * Neufeld et0al. (2010) Neufeld, J0A, Hesse, MA, Riaz, A, Hallworth, MA, Tchelepi, HA & Huppert, HE 2010 Convective dissolution of carbon dioxide in saline aquifers. Geophys. Res. Lett. 370(22).
  * Nield & Bejan (2017) Nield, D. & Bejan, A. 2017 Convection in porous media. Spinger fifth edition.
  * Novotný & Söhnel (1988) Novotný, P. & Söhnel, O. 1988 Densities of Binary Aqueous Solutions of 306 Inorganic Substances. J. Chem. Eng. Data 330(1), 49–55.
  * Ostilla-Monico et0al. (2015) Ostilla-Monico, R., Yang, Yantao, van der Poel, E.0P., Lohse, D. & Verzicco, R. 2015 A multiple-resolution strategy for Direct Numerical Simulation of scalar turbulence. J. Comput. Phys. 301, 308–321.
  * Pirozzoli et0al. (2021) Pirozzoli, S., De0Paoli, M., Zonta, F. & Soldati, A. 2021 Towards the ultimate regime in Rayleigh–Darcy convection. J. Fluid Mech. 911, R4.
  * Riaz et0al. (2006) Riaz, A, Hesse, M, Tchelepi, HA & Orr, FM 2006 Onset of convection in a gravitationally unstable diffusive boundary layer in porous media. J. Fluid Mech. 548, 87–111.
  * Sardina et0al. (2018) Sardina, G, Brandt, L, Boffetta, G & Mazzino, A 2018 Buoyancy-driven flow through a bed of solid particles produces a new form of rayleigh-taylor turbulence. Phys. Rev. Lett. 1210(22), 224501.
  * Simmons et0al. (2001) Simmons, C.0T., Fenstemaker, T.0R. & Sharp0Jr, J.0M. 2001 Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges. J. Contam. Hydrol. 520(1-4), 245–275.
  * Slim (2014) Slim, A.0C. 2014 Solutal-convection regimes in a two-dimensional porous medium. J. Fluid Mech. 741, 461–491.
  * Slim et0al. (2013) Slim, A0C, Bandi, MM, Miller, J0C & Mahadevan, L 2013 Dissolution-driven convection in a Hele–Shaw cell. Phys. Fluids (1994-present) 250(2), 024101.
  * Szulczewski et0al. (2013) Szulczewski, ML, Hesse, MA & Juanes, R 2013 Carbon dioxide dissolution in structural and stratigraphic traps. J. Fluid Mech. 736, 287–315.
  * Tsinober et0al. (2023) Tsinober, A, Shavit, U & Rosenzweig, R 2023 Numerical investigation of the influence of hydrodynamic dispersion on solutal natural convection. Water Resourc. Res. p. e2023WR034475.
  * Ulloa & Letelier (2022) Ulloa, HN & Letelier, JA 2022 Energetics and mixing of thermally driven flows in hele-shaw cells. J. Fluid Mech. 930, A16.
  * Van Der0Molen & Van0Ommen (1988) Van Der0Molen, W.H. & Van0Ommen, H.0C. 1988 Transport of solutes in soils and aquifers. J. Hydrol. 1000(1), 433 – 451.
  * van der Poel et0al. (2015) van der Poel, EP., Ostilla-Mónico, R, Donners, J & Verzicco, R 2015 A pencil distributed finite difference code for strongly turbulent wall-bounded flows. Computers & Fluids 116, 10–16.
  * Verzicco (2023) Verzicco, R 2023 Immersed Boundary Methods: Historical Perspective and Future Outlook. Annu. Rev. Fluid Mech. 550(1), 129–155.
  * Verzicco & Orlandi (1996) Verzicco, R. & Orlandi, P. 1996 A Finite-Difference Scheme for Three-Dimensional Incompressible Flows in Cylindrical Coordinates. J. Comput. Phys. 1230(2), 402–414.
  * Wang et0al. (2016) Wang, L, Nakanishi, Y, Hyodo, A & Suekane, T 2016 Three-dimensional structure of natural convection in a porous medium: Effect of dispersion on finger structure. Int. J. Greenhouse Gas Control 53, 274–283.
  * Wang et0al. (2018) Wang, L., Nakanishi, Y., Teston, A.0D. & Suekane, T. 2018 Effect of diffusing layer thickness on the density-driven natural convection of miscible fluids in porous media: Modeling of mass transport. J. Fluid Science and Technology 130(1), JFST0002–JFST0002.
  * Wang et0al. (2022) Wang, Y., Vuik, C. & Hajibeygi, H. 2022 Analysis of hydrodynamic trapping interactions during full-cycle injection and migration of co2 in deep saline aquifers. Adv. Wat. Res. 159, 104073.
  * Wen et0al. (2018) Wen, B, Akhbari, D, Zhang, L & Hesse, MA. 2018 Convective carbon dioxide dissolution in a closed porous medium at low pressure. J. Fluid Mech 854, 56–87.
  * Wettlaufer et0al. (1997) Wettlaufer, J.0S., Worster, M.0G. & Huppert, H.0E. 1997 The phase evolution of Young Sea Ice. Geophys. Res. Lett. 240(10), 1251–1254.
  * Woods (2015) Woods, AW 2015 Flow in porous rocks. Cambridge University Press.
  * Xu & Yu (2008) Xu, P & Yu, B 2008 Developing a new form of permeability and kozeny–carman constant for homogeneous porous media by means of fractal geometry. Adv. Wat. Res. 310(1), 74–81.
  * Xu et0al. (2006) Xu, X, Chen, S & Zhang, D 2006 Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv. Water Resour. 290(3), 397–407.
  * Zaman & Jalali (2010) Zaman, E & Jalali, P 2010 On hydraulic permeability of random packs of monodisperse spheres: Direct flow simulations versus correlations. Phys. A: Stat. Mech. Appl. 3890(2), 205–214.
  * Zhong et0al. (2023) Zhong, J, Liu, S & Sun, C 2023 On the thermal effect of porous material in porous media Rayleigh-Bénard convection. Flow 3, E13.
  * Zivar et0al. (2021) Zivar, D, Kumar, S & Foroozesh, J 2021 Underground hydrogen storage: A comprehensive review. Int. J. Hydrog. Energy 460(45), 23436–23462.