Convective mixing in porous media: a review of Darcy, pore-scale and Hele-Shaw studies

Abstract

Convection-driven porous media flows are common in industrial processes and in nature. The multiscale and multiphase character of these systems and the inherent nonlinear flow dynamics make convection in porous media a complex phenomenon. As a result, a combination of different complementary approaches, namely theory, simulations and experiments, have been deployed to elucidate the intricate physics of convection in porous media. In this work, we review recent findings on mixing in fluid-saturated porous media convection. We focus on the dissolution of a heavy fluid layer into a lighter one, and we consider different flow configurations. We present Darcy, pore-scale and Hele-Shaw investigations inspired by geophysical processes. While the results obtained for Darcy flows match the dissolution behaviour predicted theoretically, Hele-Shaw and pore-scale investigations reveal a different and tangled scenario in which finite-size effects play a key role. Finally, we present recent numerical and experimental developments and we highlight possible future research directions. The findings reviewed in this work will be crucial to make reliable predictions about the long-term behaviour of dissolution and mixing in engineering and natural processes, which are required to tackle societal challenges such as climate change mitigation and energy transition.

Graphic abstract

1 Introduction

A porous medium is a material consisting of a solid matrix with an interconnected void, which allows fluids to flow through it. When a fluid-saturated medium subject to the action of gravity experiences an unstable density profile, i.e. a heavy fluid parcel sitting above a less dense one, the denser fluid will eventually move and replace the lighter fluid, and vice versa. The density-driven physical mechanism inducing this motion is defined as convection, and it represents the driving force of many problems of practical interest, particularly in geophysical processes. The regular polygonally patterned crusts of salt shown in Fig. 1a, approximately a metre in diameter, are the surface signature of the vertical transport of salt, a fundamental process in arid regions. These ridges form as a result of solutal convection in the porous soil beneath the surface [1, 2]. Similarly, in supercritical geothermal systems heat supplied by a magmatic heat source produces a buoyancy-induced flow circulation due to convection [3]. Formation of sea ice (Fig. 1b) or solidification of multicomponent alloys may originate mushy layers, which consist of a porous medium filled with interstitial fluid [4]. This fluid (brine, a mixture of water and sea salt) experiences density gradients produced by differences of temperature and solute concentration, which induce convective motions within the porous layer and control the subsequent solidification dynamics [5, 6]. The above convective processes in porous media are associated with grand societal challenges, including energy transition and climate change mitigation. Understanding the underlying fluid mechanics is crucial for making reliable predictions on the evolution of the natural environment [7]. Within the many applications of importance in this context, convection in porous media has received renovated attention due to the implications it bears for geological sequestration of carbon dioxide (CO$_{2}$) [8].

Fig. 1

Examples of convection in porous media in geophysical applications. a Salt polygons at the Hoz-e Soltan (Iran) (image courtesy of [9]). These superficial formations are the result of salt-induced convective subsurface flows. b Formation of sea ice (adapted with permission from [10]). When sea ice grows, the intermediate layer between the ice exposed to the atmosphere and the ocean forms a porous solid matrix (ice) filled in the interstitial space by brine (water and salt). Salt-rich (yellow) plumes of brine drain from this mushy layer into the underlying ocean (blue). c Migration of carbon dioxide (CO$_2$) in a post-injection scenario (adapted with permission from [11]). Brine and CO$_2$ saturate the porous medium and are vertically confined by two low-permeability layers. Due to symmetry, only the right half of the reservoir is shown. (square) Dissolution of CO$_2$ in brine occurs at the interface between the currents of these fluids. (circle) Liquid phase filling the interstitial space within the pores of the rocks

Geological CO$_{2}$ storage consists of injecting large volumes of carbon dioxide in underground geological formations with the aim of permanent (or long-term) storage (Fig. 1c). These formations are typically saline aquifers and consist of a porous material confined by horizontal low-permeability layers (grey regions in Fig. 1c). The aquifers are located 1–3 km beneath the Earth surface, where the pressure is sufficient to keep the CO$_{2}$ supercritical [12, 13]. Here, a rich flow dynamics emerges. Injected CO$_{2}$ (black) is initially lighter than the fluid (brine, yellow) naturally filling the subsurface aquifer, and therefore, carbon dioxide migrates towards the upper region of the formation, driven by convection, to form a CO$_{2}$ layer that spreads horizontally. The low-permeability layer prevents CO$_{2}$ from escaping and migrating to the uppermost parts of the aquifer, from where it could eventually return to the atmosphere. At the interface between the currents of carbon dioxide and brine, the dissolution of CO$_{2}$ into the underlying brine layer takes place, originating a new mixture heavier than both starting fluids (red-to-green fluid in Fig. 1c). The dissolution process, illustrated in the squared inset of Fig. 1c, makes the interfacial layer heavier and thicker, and eventually finger-like instabilities form. The CO$_{2}$-rich solution will sink down being permanently stored in the formation. The presence of these finger-like structures makes the convective dissolution process more efficient compared to a diffusive dissolution. Such a behaviour is highly desired for CO$_{2}$ storage because injected carbon dioxide will dissolve faster preventing leakages in case of faults at the top low-permeability confining layer. In turn, the presence of nonlinear structures makes the system complex to study, and long-term predictions of the dynamics of injected carbon dioxide require huge computational efforts. An element further increasing the complexity of this scenario is represented by the finite-size pore-scale effects. At the level of the rock grains, schematically reported in the circle of Fig. 1c, the fluid moves in the interstitial space following sinuous paths, further spreading the solute transported and making predictions on the long-term behaviour even more challenging. Motivated by the CO$_{2}$ storage process, convection in porous media has been recently investigated in great detail [14]. In this work, we will review the current modelling approaches, numerical and laboratory measurements, and in particular, we will focus on the role of finite-size effects such as confinements and pore-scale dispersion.

In a convective porous medium flow, the dynamics is controlled by the relative importance of driving and dissipative mechanisms, which is quantified by the Rayleigh–Darcy number $\text {Ra}$. Convection is the driving process, and it is determined by the combination of fluid properties (density contrast), medium properties (permeability and porosity) and domain properties (gravity and domain size). Dissipative forces act against convection either as a drag force between the fluid and the solid (due to viscosity) or reducing local gradients of density (due to molecular diffusion). As a result of solute redistribution due to the tortuous fluid path in the interstitial matrix, the solute concentration field is made more uniform. This effect, labelled as dispersion, contributes as well to dissipate the potential mixing energy of the system, since the concentration gradients within the domain reduce. A key challenge in studying convective geophysical flows consists of making reliable predictions of their evolution by determining how global transport quantities, e.g. the solute flux or the mixing rate, vary as a function of $\text {Ra}$ and time. Simplified mathematical models solved numerically and theoretically have provided a clear picture of the flow behaviour at the Darcy scale [15,16,17], i.e. when a sufficiently large representative elementary volume including many pores is considered [18]. In contrast, these results disagree with the experimental measurements [19, 20], possibly suggesting that physical effects present in laboratory set-ups are not captured by the classical Darcy formulation [21].

An intuitive way to experimentally mimic a porous medium consists of filling a confined volume with solid materials, and when spherical objects are used, the medium is defined as a bead pack [22]. These experiments may be challenging, since the medium is typically hardly accessible due to its opacity, and only in recent years, non-invasive and non-intrusive measurements such as X-ray tomography and magnetic resonance imaging have become accessible [23, 24]. As a result, most of the experiments on convective flows in porous media have been performed in Hele-Shaw cells, which consist of two transparent plates separated by a narrow gap where the fluid flows [25]. The Hele-Shaw apparatus is particularly relevant because it provides optical access, and in some conditions, the flow follows a Darcy-like behaviour. In general, neither bead packs nor Hele-Shaw cells faithfully reproduce the dynamics of a Darcy flow, in which the flow structures within a porous medium are much larger than the average pore size. A difference emerged in the transport properties of bead packs experiments, Hele-Shaw experiments and numerical simulations [21]. The solute redistribution effects (dispersion) produced either by the presence of the solid obstacles in the porous matrix or by the walls in a Hele-Shaw flow have been identified as the main responsible for these discrepancies [26], and are labelled here as finite-size effects. In recent years, the advancement of theoretical, experimental and numerical techniques allowed a more precise characterisation of the flow, with accurate measurements of pore-scale dissolution rates, and a clearer picture of the influence of finite-size effects on dissolution and mixing is now available.

In this work, we review recent theoretical, numerical and experimental findings in the field of convection in porous media. This review is meant to be complementary with respect to other works [12,13,14], since we focus on dissolution and mixing with emphasis on finite-size effects. The paper is organised as follows. In Sect. 2, we describe the mathematical models and the idealised configurations used to investigate convection in porous media, and we derive a unified formulation to evaluate and relate mixing in different flow configurations. In Sects. 3 and 4, we review the results obtained in Rayleigh–Bénard and one-sided configurations, respectively. Finite-size effects possibly leading to the discrepancy observed between experiments and simulations are discussed in Sect. 5. Finally, in Sect. 6 we summarise the results discussed and present recent experimental developments, together with an overview of additional effects not present in the configurations discussed in Sects. 3 and 4.

2 Modelling of convection

2.1 Pore-scale modelling

Convective flows are produced by the presence of unstable density gradients within an accelerated fluid domain. These density differences drive the flow towards a more stable configuration, decreasing the gravitational potential energy within the system [27]. We consider problems in which convection is induced by the presence of a scalar quantity (e.g. solute concentration or temperature) that modifies the density field of the flow. For simplicity, in this review we will define the parameters in case of solute convection, but the findings extend to the case of thermally driven convection unless explicitly mentioned.

The maximum density difference within the domain, $\Delta \rho$, determines the strength of the convective flow. On the other hand, (molecular or thermal) diffusion reduces the local scalar gradients diminishing the driving force of the flow, and viscosity is responsible for energy dissipation due to friction. In a free fluid (i.e. in the absence of a porous medium), the relative importance of these two contributions is quantified by the Rayleigh number $\text {Ra}_{T}$ defined on the characteristic length scale of the flow H

$$\begin{aligned} \!\text {Ra}_{T}=\frac{g\Delta \rho H^3}{\mu D} , \end{aligned}$$

(1)

where g is the acceleration due to gravity, D is the molecular diffusivity and $\mu$ is the fluid dynamic viscosity. The ratio of kinematic viscosity to solute diffusivity D (or molecular diffusivity) determines the Schmidt number

$$\begin{aligned} \text {Sc}=\frac{\mu }{\rho _r D}, \end{aligned}$$

(2)

with $\rho _r$ the average (or reference) fluid density within the domain. Similarly, for thermally driven flows one can define the Prandtl number ($\text {Pr}$), in which the molecular diffusion is replaced by its thermal counterpart.

Modelling heat or mass transport at the pore scale requires to resolve the flow within the interstitial space. Momentum transport is controlled by continuity and Navier–Stokes equations, respectively:

$$\begin{aligned}{} & {} \nabla \cdot {\tilde{\textbf{u}}}=0\text {,} \end{aligned}$$

(3)

$$\begin{aligned}{} & {} \tilde{\rho }\left[ \frac{\partial \tilde{\textbf{u}}}{\partial t}+\left( {\tilde{\textbf{u}}}\cdot \nabla \right) \tilde{\textbf{u}}\right] =-\nabla \tilde{p} +\mu \nabla ^2{\tilde{\textbf{u}}} + \tilde{\rho }\textbf{g}\text {,} \end{aligned}$$

(4)

where $\tilde{\textbf{u}}$, ${\tilde{\varvec{\rho }}}$ and $\tilde{p}$ are the velocity, density and pressure fields, respectively, and $\textbf{g}$ indicates acceleration due to gravity. We assumed that the Boussinesq approximation applies, which is reasonable for geophysical processes such as carbon sequestration [28]. (Additional details on this assumption are provided in Sect. 6.2.) The fluid density $\tilde{\rho }$ is typically defined by an equation of state (EOS) that depends on both solute concentration and fluid temperature. (Other scalars present in the system may be similarly treated.) When linearised, the EOS may be rewritten to obtain the density $\tilde{\rho }$ as a function of temperature ($\tilde{T}$) and concentration ($\tilde{C}$) (possible limitations of this approach are discussed in Sect. 6.2). With respect to the value $\tilde{\rho }_r$ defined at the reference state ($\tilde{C}_r,\tilde{T}_r$), it reads

$$\begin{aligned} \tilde{\rho }=\tilde{\rho }_r+\alpha _c(\tilde{C}-\tilde{C}_r)+\alpha _t(\tilde{T}-\tilde{T}_r), \end{aligned}$$

(5)

where $\alpha _c,\alpha _t$ are the expansion coefficients relating the density to the variations of concentration and temperature, respectively, being typically $\alpha _c>0$ and $\alpha _t<0$. In this case, assuming the presence of solute scalars only, Eq. (5) reduces to the form $\tilde{\rho }=\tilde{\rho }(\tilde{C})=\tilde{\rho }_r+\alpha _c(\tilde{C}-\tilde{C}_r)$.

Solute conservation is accounted by the advection–diffusion equation:

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

(6)

Equations (3)–(6) are solved for the fluid domain to determine the evolution of the flow at the pore scale [29]. When heat transport is considered, a diffusive heat flux in the solid matrix may be also accounted [30, 31] (additional details will be provided in Sect. 3.2). The presence of additional phases is not discussed in this review, and we refer to [32] for pore-scale modelling approaches of multiphase flows.

Notwithstanding recent significant improvements of numerical schemes and computational infrastructures, resolving real-scale convective flows from the pore level to the reservoir scale requires a computational effort that is beyond the present capabilities. To overcome this issue, a possible approach consists of modelling the flow at an intermediate scale between pores length scale and domain height, i.e. the Darcy scale. Despite missing a precise description of the flow dynamics at the pore level, Darcy models have proved to be a reliable framework to determine the overall long-time behaviour for the transport of species in convective porous media flows [2, 10, 12, 33]. In the following, we will describe under which assumptions a convective flow in porous media can be modelled as a continuum via the Darcy flow approximation.

2.2 Darcy model and dispersion

A possible strategy to model flows in porous media consists of taking the average of relevant quantities (velocity, concentration and pressure fields) over a representative volume that contains several pores [18]. An illustrative example is sketched in Fig. 2. The size of the volume (indicated as representative elementary volume, REV) over which the average is computed is larger compared to the pores length scale d, but still smaller than the domain reference length H. The Darcy model is based on empirical observations initially proposed more than 150 years ago [34], and the later derived analytically by Whitaker [35]. We refer to [18, 36] and references therein for additional details. The key assumption of the Darcy equation is that the average flow velocity over the representative volume is proportional to the pressure gradient applied to the volume via the fluid viscosity and a property of the medium defined as permeability. These conditions are achieved when the flow inertia is negligible compared to viscous forces [18]. In the following, we will characterise the medium properties and the governing flow parameters, and finally, we will discuss a model for the Darcy flow.

Fig. 2

Model of the flow at different scales. a At the Darcy level, all flow quantities are obtained as averaged over the REV. Solid boundaries (in this example at the bottom of the domain) are impermeable to fluid, i.e. the velocity component perpendicular to the wall is zero ($\textbf{n}$ is the unit vector normal to the wall). However, slip along this boundary is possible. b At the pore level, the fluid phase flows within the interstitial space of the solid matrix, which is made of impermeable solid objects. Over the surface of each of these objects, no-slip boundary condition applies

The characteristic geometrical properties of the solid matrix and its intimate interaction with the interstitial fluid determine the flow behaviour. The main macroscopic parameters used to characterise a porous medium are: (i) porosity $\phi$, defined as the ratio of volume of fluid to the total volume (fluid and solid), and (ii) permeability k, a measure of the resistance opposed by the medium to the flow. For a given porous medium, the Darcy number represents the relative importance of permeability and its cross-sectional reference area. With respect to the domain reference length scale H, the Darcy number reads:

$$\begin{aligned} \text {Da}=\frac{k}{H^2}. \end{aligned}$$

(7)

A convective flow is driven by density differences, and therefore a possible velocity scale is the buoyancy velocity U, i.e. the free fall velocity of a parcel of immiscible fluid surrounded by fluid having a density contrast $\Delta \rho$, which is defined as

$$\begin{aligned} U=\frac{g\Delta \rho k}{\mu } . \end{aligned}$$

(8)

We observe that U is independent of the any length scale, and it relates to the fluid ($\Delta \rho ,\mu$), medium (k) and domain (g) properties. In addition to the domain length scale (H), one can consider as a reference length scale the distance $\ell$ over which advection and diffusion balance [17]

$$\begin{aligned} \ell =\frac{\phi D}{U} \end{aligned}$$

(9)

(in thermal convection, $\ell =D/U$ with D representing the thermal diffusivity). The evolution of the fluid layer is controlled by buoyancy, which tends to drive the flow towards a stable configuration, and diffusion, acting to reduce local concentration gradients and increasing the mixing of solute in the domain. The relative importance of the strength of these contributions is evaluated by the Rayleigh–Darcy number $\text {Ra}$

$$\begin{aligned} \!\text {Ra}=\frac{H}{\ell }\text {,} \end{aligned}$$

(10)

obtained combining the Rayleigh number $\text {Ra}_{T}$ (1) and the Darcy number $\text {Da}$ (7). In particular, $\text {Ra}=\!\text {Ra}_{T}\text {Da}/\phi$ in the instance of solutal convection, with the solid being impermeable to the solute fluxes, and $\text {Ra}=\!\text {Ra}_{T}\text {Da}$ for thermally driven cases in conductive media. While in the thermal case an equilibrium between the solid and the fluid phases may be achieved, in solutal convection the solid phase is always solute-free. Notwithstanding this difference, when thermal equilibrium locally occurs between the solid and the fluid phases, results for thermal convection can be equally interpreted as results for solutal convection, and vice versa, provided that the Rayleigh–Darcy number is matched [14]. The Rayleigh–Darcy number includes all the macroscopic properties of the system: domain (g, H), medium ($k,\phi$) and fluid ($D,\mu ,\Delta \rho$) properties. In addition, when the spatial coordinates are made dimensionless with respect to $\ell$ (9), $\text {Ra}$ can be interpreted as the dimensionless domain height [17].

A Darcy-type flow occurs when the size of the flow structures is much greater than the reference length of the REV [14]. The reference length scale is in this case the pore scale, which is proportional to $\sqrt{k}$. In quantitative terms, the criterion above is fulfilled when: (i) the pore-scale Reynolds number is small, i.e. viscous dissipation $(\mu U)$ dominates over inertia $(\rho _{r}U^{2}\sqrt{k})$, and (ii) the smallest length scale of the flow ($\ell$) is large compared to the pore size ($\sqrt{k}$). These conditions are matched if:

$$\begin{aligned} \frac{\rho _{r}U^{2}\sqrt{k}}{\mu U}\ll 1\Rightarrow \text {Re}&=\frac{\rho _{r}U\sqrt{k}}{\mu }=\frac{\!\text {Ra}\text {Da}^{1/2}}{\textrm{Sc}}\ll 1 \end{aligned}$$

(11)

$$\begin{aligned} \frac{\ell }{\sqrt{k}}\gg 1\Rightarrow \text {Pe}&=\frac{\sqrt{k}}{\ell }=\!\text {Ra}\text {Da}^{1/2}\ll 1, \end{aligned}$$

(12)

i.e. when Reynolds ($\text {Re}$) and Péclet ($\text {Pe}$) numbers are much less than unity. Note that in these definitions the pore length scale ($\sqrt{k}$) and the buoyancy velocity (U) are used as length and velocity scales, respectively.

We consider a fluid-saturated homogeneous and isotropic porous medium with porosity $\phi$ and permeability k (Fig. 2a) fulfilling the conditions (11)–(12). The flow field is fully described by the continuity and Darcy equations, respectively:

$$\begin{aligned}{} & {} \nabla \cdot \textbf{u}=0\quad \end{aligned}$$

(13)

$$\begin{aligned}{} & {} \textbf{u}=-\frac{k}{\mu }\left( \nabla p+\rho \textbf{g}\right) \text {.} \end{aligned}$$

(14)

Note that in this case $\textbf{u}$ is the seepage or Darcy velocity, and it represents the value of fluid velocity averaged over the REV (Fig. 2b). It is related to the fluid velocity averaged over the fluid phase of the REV ($\tilde{\textbf{u}}$) via the Dupuit–Forchheimer relationship $\textbf{u}=\phi \tilde{\textbf{u}}$ [18]. Same applies to pressure p and density $\rho$.

The evolution of the concentration field is controlled by the advection–diffusion equation

$$\begin{aligned} \phi \frac{\partial C}{\partial t}+\nabla \cdot (\textbf{u}C-\phi D \nabla C)=0 , \end{aligned}$$

(15)

where t is time, C is the concentration averaged over the REV and D is the solute diffusivity, which is assumed constant and independent from the flow. In a more general formulation discussed in Sect. 2.2.1, this coefficient may be replaced by a dispersion tensor $\textbf{D}$ that depends on the local flow conditions ($\textbf{u}$) or the fluid properties (Sc). While the solid is commonly impermeable to solute, in the thermal case a diffusive heat flux may occur within the solid matrix. In case of thermal equilibrium between the solid and the liquid phases, Eqs. (13)–(15) keep being valid.

2.2.1 Dispersion

Solute redistribution induced by the fluid carrying the solute and flowing through the porous medium is defined as dispersion [22]. This mechanism, which has the effect of homogenising the solute concentration field, adds to the contribution of molecular diffusion. For this reason, these two contributions that are originated by very different physical mechanisms are often grouped within a unique formulation. In porous media, dispersion may arise due to several reasons: pore-scale change of flow direction (mechanical dispersion), heterogeneous permeability fields (large-scale dispersion) or other mechanisms, such as no-slip at the boundary of the pores or dead-end pores (anomalous dispersion). These effects are the result of the pore-scale dynamics and can be considerably more effective (up to few orders of magnitude) than the solute spreading due to molecular diffusion. Therefore, it may be necessary to account for the presence of dispersion when modelling the flow at the Darcy scale. Here we consider the contribution of mechanical dispersion and molecular diffusion, usually grouped in a term defined as hydrodynamics dispersion. For simplicity, in the following we will indicate this mechanism as dispersion, ad we refer to [37] for a general theoretical discussion on dispersion-induced mixing.

A classical approach to account for the effects of dispersion consists of replacing the molecular diffusion coefficient, D in Eq. (15), with a dispersion tensor $\textbf{D}$ which depends on the local flow conditions. Typically, the dispersion tensor is anisotropic and aligned with the flow, meaning that it can be decomposed into two components in the directions parallel ($D_{\textrm{L}}$, longitudinal dispersion) and perpendicular ($D_{\textrm{T}}$, transverse dispersion) to the local flow velocity $\textbf{u}$. This model is labelled as Fickian dispersion [38]. With these assumptions, the dispersion tensor takes the form:

$$\begin{aligned} \textbf{D} = D\textbf{I}+(\alpha _{\textrm{L}}-\alpha _{\textrm{T}})\frac{\textbf{u}\textbf{u}}{\vert \textbf{u}|} +\alpha _{\textrm{T}}\textbf{u}\textbf{I}, \end{aligned}$$

(16)

where $\textbf{I}$ is the identity tensor, and the coefficients $\alpha _{\textrm{L}}=D_{\textrm{L}}/U$ and $\alpha _{\textrm{T}}=D_{\textrm{T}}/U$ correspond to the dispersivities of the medium in longitudinal and transverse directions, respectively. For solute transport and $\text {Pe}\gg 1$, dispersion in the cross-flow direction is typically 1 order of magnitude smaller than in the stream flow direction (possible limitations of this assumption are discussed in [39]). The ratio of these two contributions is quantified by the dispersivity ratio

$$\begin{aligned} r=\frac{D_\textrm{L}}{D_\textrm{T}}. \end{aligned}$$

(17)

The magnitude of $D_\textrm{L}$ and $D_\textrm{T}$ is estimated with the aid of correlations based on experiments and simulations (see [22], and references therein). Longitudinal and transverse dispersion coefficients depend on many parameters, namely Schmidt number [40], Reynolds number [41], tortuosity of the medium [42], Péclet number [43] and fluid phases [42]. We refer to [44] for a review of numerical, experimental and theoretical works in this field.

We consider here an example of a medium composed of monodispersed beads, typical of numerical and experimental set-ups commonly employed in pore-scale investigations, and we show that the dispersivity ratio r may considerably vary as a function of the Péclet number of the flow (a similar procedure applies for different media—porosity, tortuosity—and flow—Péclet number—properties). We consider the case of solutal convection in a monodispersed bead pack at $\text {Sc}\le 550$. For a monodispersed close random packing, the porosity is $\phi = 0.37$ [45, 46] and the tortuosity (the ratio of actual flow path length to the straight distance between the ends of the flow path [36]) is $\tau =0.68$ (see [47], and references therein). We use the empirical correlations proposed by Delgado [42], obtained for laboratory experiments, i.e. in architecture-controlled media, whereas we refer to [48] for a review of dispersion relations in field-scale data. The results proposed by Delgado [42] are valid for liquids and at $\text {Sc}\le 550$, and we report the dispersivity ratio in Fig. 3 (for $\text {Sc}>550$, similar correlations are provided). Four main flow regimes have been identified, for increasing Pe: (i) diffusion regime, with molecular diffusion being the dominant mechanism; (ii) diffusion and mechanical dispersion, when the two contributions are comparable; (iii) pure mechanical dispersion, when the influence of molecular diffusion is negligible; and (iv) non-Darcy, when the effects of inertia and turbulence cannot be neglected. Note that these correlations are obtained from experimental measurements performed for a wide parameters space, and a sharp separation between these regimes is hard to identify. A theoretical prediction is available for low ($r=1$, [49]) and high ($r=6$, [42] and references therein) Péclet numbers. A similar regime classification has been also proposed by Perkins and Johnston [50].

With this example, we have shown that in general r varies with Pe and Sc among the other parameters, and also across the scales [51]. To simplify the picture, a possible approach used in numerical simulations consists of fixing the values of $D_{\textrm{L}}$ and $D_{\textrm{T}}$ or r [26], which is a reasonable approximation if a narrow range of Pe is considered. Results on the effect of dispersion on convective flow are presented in Sect. 5.3.

Fig. 3

Dispersivity ratio $r=D_\textrm{L}/D_\textrm{T}$ shown for porosity $\phi =0.37$, tortuosity $\tau =0.68$ and different Schmidt numbers, namely $\text {Sc}= 50$, 150, 250, 350, 450 and 550. The correlations proposed by Delgado [42] have been employed. The advective flow is divided in several regimes, discussed in the text

2.3 Flow configurations and quantification of mixing

Convective processes of practical interest are characterised by the mixing of one or more scalar quantities (e.g. the concentration of a dispersed solute phase) in the ambient fluid and predicting the time required to achieve a certain degree of mixing may be necessary. In the instance of geological carbon sequestration, for example, it is desired to find the time required to dissolve a considerable fraction of the CO$_{2}$ injected, to assess the reliability of a given sequestration site. These estimates can be obtained via experiments and simulations in representative flow configurations, which are well controlled and designed to reproduce the main features observed in environmental and industrial cases. In this section, we will first introduce the main flow configurations investigated in literature, with clear indication of the initial and boundary conditions. Then we will define a general framework and identify relevant observables required to quantify the mixing and analyse the evolution of the system.

Fig. 4

Flow configurations (adapted with permission from [52]). a Sketch of boundary conditions applied at the top (label 1) and bottom (label 0) boundaries in terms of flux F and concentration C. All boundaries are impermeable to fluid ($\textbf{u}\cdot \textbf{n}=0$), and side boundaries may be also considered as periodic. The reference frame (x, z) and gravity ($\textbf{g}$) are also indicated. Three flow configurations are shown: b Rayleigh–Bénard, c one-sided and d Rayleigh–Taylor. An exemplar field obtained for two-dimensional simulations at $\text {Ra}=7244$ is reported. The field is taken at the time indicated by the green arrows in panels (b-ii), (c-ii) and (d-ii), where the evolution of the parameters $\widehat{\chi }$, $\widehat{F}$ and $d_{t}\langle C^{2}\rangle /2$ is reported (the operator $d_{t}$ stands for the time derivative). Quantities are computed as in Eq. (25) and made dimensionless with respect to the length scale $\mathcal {L}=\ell$. The time-averaged value of $\widehat{F}$ is also shown (dashed lines) in panels (b-ii) and (c-ii)

Three archetypal flow configurations are generally employed to investigate the dynamics of convection in porous media. They consist of analogue systems that help us to have a comprehension of specific scenarios occurring in nature. A sketch to illustrate possible boundary conditions applied is shown in Fig. 4a. At the top (label 1) and bottom (label 0) boundaries, both flux F and concentration C may be prescribed. (The flux will be defined more precisely later in this section.) All boundaries are considered impermeable to fluid, i.e. no penetration condition applies ($\textbf{u}\cdot \textbf{n}=0$, being $\textbf{u}$ the fluid velocity and $\textbf{n}$ the vector perpendicular to the boundary. However, periodic conditions on the side boundaries may be also considered for convenience in numerical studies, with no difference on the modelling described in the following. In all cases considered here, the domain boundaries are assumed impermeable to fluid, and the fluid is supposed initially still [$\textbf{u}(t=0)=0$]. The boundary conditions for the solute (fixed concentration or no flux) will determine the nature of the system considered (steady or transient), whereas the initial condition for the solute (uniform concentration, or two fluid layers with different concentration) will control the flow evolution. The flow configurations considered are:

• i. Rayleigh–Bénard (Fig. 4b-i): the solute concentration is fixed at the horizontal boundaries, so that the density of the fluid at the bottom wall ($C=C_{0}$) is lighter than the density of the fluid at the top wall ($C=C_{1}$) [53, 54]. This unstable flow attains a statistically steady state, which is rigorously steady for sufficiently low Rayleigh–Darcy numbers [55]. A scalar flux is possible through the upper ($F=F_{1}$) and the lower ($F=F_{0}$) boundaries.

• ii. One-sided (Fig. 4c-i): the concentration is imposed at the upper wall, where a solute flux is also possible ($C=C_{1}, F=F_{1}$), and the domain is impermeable to solute at the lower wall ($\partial C/\partial z=0$, corresponding to $F_{0}=0$). This configuration originates a time-dependent flow, and the domain, initially filled with uniform solute concentration $C=C_0$, is gradually filled with the solute coming from the upper boundary [15, 17, 56].

• iii. Rayleigh–Taylor (Fig. 4d-i): both walls are impermeable to the scalar ($F_{0}=F_{1}=0$). The domain initially consists of two uniform layers of different density ($C=C_{1}$ for the upper portion, and $C=C_{0}$ for the lower portion), so that the flow configuration is unstable [57, 58]. Solute mixing evolves controlled by the dynamics of the flow structures.

The flow configurations considered differ in terms of boundary conditions and evolution, and suitable flow observables are required to estimate the mixing state of each system. For instance, the Sherwood number $\text {Sh}$, defined as the ratio of the convective to the diffusive mass transport, is suitable in solute-permeable domains (e.g. the Rayleigh–Bénard case), but it does not provide any indication in closed domains (e.g. the Rayleigh–Taylor case). Therefore, in each flow configuration different quantities are used, which are related through exact mathematical relations that are derived here. Following [21], we take the advection–diffusion equation (15) multiplied by C, and we integrate over the entire domain. We use the hypothesis of incompressibility of the flow (14) together with the impermeability of the boundaries to the fluid. (Note that the same result is achieved assuming periodicity in horizontal direction.) After some algebraic manipulations, we obtain the following exact global relation:

$$\begin{aligned} \frac{\phi }{2}\frac{\textrm{d} \langle C^{2}\rangle }{\textrm{d} t} = \frac{\phi }{H}\left( C_{1}F_{1}+C_{0}F_{0}\right) -\phi \chi , \end{aligned}$$

(18)

where $\langle \cdot \rangle$ indicates the volume average. Equation (18) relates the mean squared concentration, the solute flux through the walls F and the mean scalar dissipation within the domain $\chi$, respectively, defined as

$$\begin{aligned} F_i=\frac{D}{L}\int _{0}^{L}\frac{\partial C}{\partial z}\biggr |_{z=z_i}\text {d}x\quad {\text {with}}\, i=\{0,1\} , \end{aligned}$$

(19)

with L domain width, and

$$\begin{aligned} \chi = D\langle |\nabla C|^{2}\rangle . \end{aligned}$$

(20)

When C is defined as a mass concentration, F may be interpreted as the average mass of solute that enters (or leaves) the domain per unit of surface area and time. Equation (18) can be interpreted as follows. The rate of change of mean squared concentration within the domain is the result of external contributions ($F_{0},F_{1}$, either positive or negative) and dissipation of mixing energy. ($\chi$ is always positive; therefore, it contributes to a reduction of scalar variance $\langle C^{2}\rangle$.)

To enable comparisons among different systems, a possible set of dimensionless variables consists of $\mathcal {L}$ for lengths, $\phi \mathcal {L}/U$ for time, and U for velocities. The concentration C is made dimensionless as $\widehat{C}=(C-C_0)/\Delta C$, where $\widehat{\cdot }$ indicates dimensionless quantities and $\Delta C = C_1-C_0$. Making Eq. (15) dimensionless with these variables and proceeding as above, we obtain a dimensionless form of Eq. (18) that reads:

$$\begin{aligned} \frac{1}{2}\frac{\textrm{d} \langle \widehat{C}^{2}\rangle }{\textrm{d} \widehat{t}} = \frac{\phi D}{U\mathcal {L}}\left( \widehat{F}-\widehat{\chi }\right) \end{aligned}$$

(21)

with $\widehat{F}=F_1\mathcal {L}/(D\Delta C)$ the dimensionless flux and $\widehat{\chi }=\chi \mathcal {L}^2/[D(\Delta C)^{2}]$ the dimensionless mean scalar dissipation. Note that in this expression the contribution of the flux at the bottom boundary vanishes, due to the set of dimensionless variables considered. The reference length scale $\mathcal {L}$ has not been defined yet and it can be conveniently set in each configuration. With respect to the systems previously introduced, the following scenarios appear:

• i. Rayleigh–Bénard (Fig. 4b-ii): after an initial transient phase, the system attains a statistically steady state [54, 59]. The time average of Eq. (21) returns

  $$\begin{aligned} \overline{\widehat{F}} = \overline{\widehat{\chi }}, \end{aligned}$$

  (22)

  where $\overline{\cdot }$ indicates the time-averaging operator. We observe in Fig. 4b-ii that while a nonzero instantaneous contribution $\textrm{d} \langle \widehat{C}^{2}\rangle /\textrm{d} \widehat{t}$ is present, $\overline{\widehat{F}}$ and $\overline{\chi }$ fluctuate around their time-averaged value (black dashed line). Note that the reference length scale normally used in this configuration is $\mathcal {L}=H$, which gives in (21) the prefactor $\phi D/(U\mathcal {L})=1/\!\text {Ra}$. The quantity used to evaluate the mass transfer in this configurations is the Sherwood number

  $$\begin{aligned} \text {Sh}=\frac{H}{\Delta C L}\overline{\int _{0}^{L}\frac{\partial C}{\partial z}\biggr |_{z=z_1}\text {d}x} , \end{aligned}$$

  (23)

  defined as the relative contribution of convective and diffusive to diffusive mass transport. Using the definition of $\widehat{F}$ and Eq. (22), $\text {Sh}$ can be related to the flux and the mean scalar dissipation [53]:

  $$\begin{aligned} \text {Sh}=\!\text {Ra}\overline{\widehat{F}}=\!\text {Ra}\overline{\chi }. \end{aligned}$$

  (24)

• ii. One-sided (Fig. 4c-ii): the domain is impermeable to solute at the lower wall ($F_{0}=0$). By setting $\mathcal {L}=\ell$ as defined in (9), Eq. (21) is independent of $\text {Ra}$ and reads

  $$\begin{aligned} \frac{1}{2}\frac{\textrm{d} \langle \widehat{C}^{2}\rangle }{\textrm{d} \widehat{t}} = \widehat{F}-\widehat{\chi }, \end{aligned}$$

  (25)

  where

  $$\begin{aligned} \widehat{F}=\frac{\phi D}{U\Delta C}\frac{1}{L}\int _{0}^{L}\frac{\partial C}{\partial z}\biggr |_{z=z_1}\text {d}x . \end{aligned}$$

  (26)

  This choice for $\mathcal {L}$ is convenient to compare systems having different $\text {Ra}$ because the value of $\widehat{F}$ appears to be universal, as will be later discussed in Sect. 4. The time-dependent flow originated from this configuration consists of three main flow regimes [17, 60, 61]. Initially ($\widehat{t}<10^{3}$) diffusion dominates and a high-concentration high-density unstable fluid layer thickens. At a later stage ($\widehat{t}<16\!\text {Ra}$), convection takes place and plumes formed at the top boundary layer grow and invade the domain. In this phase $\widehat{F}$ is statistically steady and characterised by a value (black dashed line) that is independent of the Rayleigh–Darcy number considered. A similar behaviour holds for $\chi$, but a closer inspection reveals that after the fingers reach the bottom ($\widehat{t}>10\!\text {Ra}$) an increase of $d \langle \widehat{C}^{2}\rangle /d \widehat{t}$ is observed. A corresponding decreasing behaviour is reflected in $\chi$, but with half the amplitude. After the upper layer of the domain is also saturated ($\widehat{t}>16\!\text {Ra}$), the dissolution flux $\widehat{F}$ drops, and the system enters the shutdown regime.

• iii. Rayleigh–Taylor (Fig. 4d-ii): the domain is impermeable to the solute ($\widehat{F}=0$) and Eq. (21) reads

  $$\begin{aligned} \frac{1}{2}\frac{\textrm{d} \langle \widehat{C}^{2}\rangle }{\textrm{d} \widehat{t}} =-\frac{\phi D}{U\mathcal {L}} \widehat{\chi }. \end{aligned}$$

  (27)

  The flow is initialised considering two fluid layers of different density in an unstable configuration. Equation (27) suggests that all the potential energy initially stored by keeping the two phases segregated is dissipated as time evolves. Both $\ell$ and H can be considered as reference length scales, depending on which part of the flow evolution is considered. However, De Paoli et al. [52] have shown that $\mathcal {L}=\ell$ provides a universal picture for the evolution of $\widehat{\chi }$, and results are presented in Fig. 4d-ii using this length scale. Similarly to what observed in the one-sided configuration, the flow is initially controlled by diffusion ($\widehat{t}<10^{3}$). Afterwards ($10^{3}<\widehat{t}<\!\text {Ra}/2$) the formation of fingers is observed, which merge and grow, accelerating mixing. In this phase, occurring at $\text {Ra}<\widehat{t}<3\ \!\text {Ra}$ in the simulations considered, $\widehat{\chi }$ is observed to increase in Darcy simulations, as shown in Fig. 4d-ii, whereas it decreases in pore-scale simulations, due to finite-size effects [29]. The limits in which these regimes set in are indicative, as the flow evolution in strongly influenced by the initial perturbation. When the domain is nearly saturated, a stable density profile is achieved, and local concentration gradients are not sufficient to sustain convection, which is in turn overcome by diffusion. Correspondingly, scalar dissipation is observed to reduce, asymptotically attaining a zero value in correspondence of a uniformly mixed domain.

A major proportion of recent studies focused on the determination of correlations of the mixing parameters ($\widehat{F}, \text {Sh}$ or $\widehat{\chi }$) with the flow parameter ($\text {Ra}$). These results will be reviewed in Sects. 3 and 4 for the Rayleigh–Bénard and the one-sided configurations, respectively.

3 Rayleigh–Bénard convection

Rayleigh–Bénard convection produces the statistically steady flow discussed in Sect. 2.3, with mass transfer properties quantified by the Sherwood number, a time-averaged ratio of total (convective and diffusive) to diffusive mass transport at the boundaries of the domain defined in Eq. (23). Alternatively, the Nusselt number is used in case of thermal convection. In this section, we will review the results relative to Darcy and pore-scale flows in this configuration.

3.1 Darcy flow

In the Darcy case [Eqs. (13)–(15)], the system is uniquely controlled by the Rayleigh–Darcy number $\text {Ra}$, defined in Eq. (10), which sets the flow structure. The behaviour of $\text {Sh}$ with $\text {Ra}$ is reported in Fig. 5 for Darcy studies available in literature for two- [31, 54, 62, 63] and three-dimensional [59, 64] simulations. We briefly recall here the main features of the flow, and we refer to [14] for a detailed review of the flow structure. For $\text {Ra}<4\pi ^{2}$, the mass transport is purely diffusive [65, 66] and no convective motion arises ($\text {Sh}=1$). The flow is maintained quiescent by the dissipative (diffusive) effects that dominate over convection. For increasing $\text {Ra}$, instabilities appear in the form of steady convective rolls [55] with corresponding increase of the convective mass transfer. When $\text {Ra}\approx 400$, unsteady boundary layer instabilities take place and become progressively dominant. When the driving force is sufficiently large, namely at $\text {Ra}\approx 1300$ and $\text {Ra}\approx 1700$ for two- and three-dimensional systems [53, 64], respectively, these instabilities turn into a dynamic formation of small plumes at the boundary layer, which eventually grow and merge into larger plumes spanning the entire domain height. In this stage, the flow enters the high-$\text {Ra}$ regime [54]. The dynamics described above is similar in two- and three-dimensional domains. However, in the three-dimensional case the flow pattern obtained at low $\text {Ra}$ may be affected by the initial condition, i.e. different flow structures are obtained starting from different initial concentration distributions. In addition, hysteresis effects have been observed in the two-dimensional case [53]: when the flow is initialised using a solution obtained at higher $\text {Ra}$, the flow structure (number of rolls) and the transport properties ($\text {Sh}$) differ from those obtained starting from, for example, a linear temperature distribution or from a solution obtained at lower $\text {Ra}$. As a starting point, $\text {Ra}=1255$ is used by Otero et al. [53], and $\text {Ra}$ is progressively decreased. The system evolves following two distinct branches (Fig. 5a), both differing from the solution obtained for increasing $\text {Ra}$.

Fig. 5

Darcy simulations. panel a Sherwood number ($\text {Sh}$) as a function of Rayleigh–Darcy number ($\text {Ra}$), and panel b in compensated form ($\text {Sh}/\!\text {Ra}$). Results reported are obtained in two- [31, 54, 62, 63] and three-dimensional [59, 64] simulations of homogeneous and isotropic porous media. Best fitting laws at high Rayleigh–Darcy numbers for two- (black solid line) and three-dimensional flows (red solid line), respectively Eqs. (28) and (29), are also reported. In panel a a subset of the two-dimensional data of [53] (blue stars) is also shown to mark the presence of hysteresis effects. The solution obtained at $\text {Ra}=1255$ was used as initial condition for these simulations. As $\text {Ra}$ is decreased, the system evolves on two branches (blue lines), both differing from the solution obtained for increasing $\text {Ra}$

Determining the scaling of $\text {Sh}$ with $\text {Ra}$ in the high-$\text {Ra}$ regime has been object of active investigation in recent years, also due to the improvements of computational capabilities. In the frame of free fluids (i.e. no porous medium), Malkus [67] and Howard [68] proposed that at sufficiently high Rayleigh numbers the interior of the domain is well mixed, and the temperature gradients are localised at the wall boundary layers. The Sherwood number is then obtained as a result of the diffusive heat flux across these layers, which is inversely proportional to their thickness, and for porous media, it is predicted to scale linearly with $\text {Ra}$, $\text {Sh}\sim \!\text {Ra}$. An accurate phenomenological description of the flow and scaling arguments is provided by Hewitt [14]. The linear scaling best fitting the two-dimensional numerical results [54] is

$$\begin{aligned} \text {Sh}= 0.0069 \!\text {Ra}+2.75, \end{aligned}$$

(28)

and it agrees also with the best known theoretical upper bound, for which $\text {Sh}\le 0.0297\!\text {Ra}$ [53]. The asymptotic scaling proposed by Hewitt et al. [54] [solid black line in Fig. 5a] fits well the numerical results, and it is in agreement with the above mentioned linear predictions: The compensated Sherwood number (Fig. 5b) approaches in this case the asymptotic value (0.0069).

In three-dimensional domains, the situation differs, as the compensated Sherwood number has not reached yet the asymptotic linear scaling (Fig. 5b). The best fitting is in this case provided by Pirozzoli et al. [59]

$$\begin{aligned} \text {Sh}=0.0081\!\text {Ra}+ 0.067\!\text {Ra}^{0.61}, \end{aligned}$$

(29)

which consists of a linear relation with sublinear corrections. The discrepancy existing between the scaling obtained in three-dimensional porous media and the linear asymptotic prediction for $\text {Ra}\rightarrow \infty$ is due to the different flow structure produced by the additional degree of freedom provided by the third spatial dimension, i.e. the flow has not reached yet the asymptotic state. It was estimated [59] that in three-dimensional domains the asymptotic regime sets in at $\text {Ra}\approx 5 \times 10^5$, i.e. more than one order of magnitude beyond the threshold identified in two-dimensional flows, and further investigations at $\text {Ra}\approx 10^{6}$ are required to confirm this finding.

Resolving the flow equations at the Darcy scale at large $\text {Ra}$ may require extensive computational resources [58, 59]. An interesting approach proposed to overcome this obstacle consists of a new modelling strategy labelled as large-mode simulation (LMS) [69]. With the aid of a scale analysis, Jenny et al. [69] observed that: (i) large-scale structures are responsible for the bulk of the production of concentration variance, (ii) variance dissipation is dominated by the small diffusive scales, and (iii) both production and dissipation rates are independent of the Rayleigh–Darcy number. On this ground, they propose a LMS model in which closure is achieved by replacing the actual diffusivity with an effective one, in analogy with large eddy simulations for turbulent flows. LMS is based on resolving the low-wave number dynamics only, whereas the effect of the unresolved scales on the large ones is modelled. Results obtained with this new strategy are promising to enable simulations for long-term predictions of convective porous media flows in practical settings.

3.2 Pore-scale flow

Recent developments in computational methods allowed numerical solution of pore-resolved convective flow models, defined by Eqs. (3)–(6). Unlike the Darcy case, in pore-scale problems the flow properties cannot be lumped into a single governing parameter, and the contribution of several flow features has to be considered. With respect to the medium, obstacles shape and arrangement determine the medium permeability. The medium conductivity influences heat transport through the solid phase, and the volume fraction of solid sets the porosity. Concerning the fluid and the scalar transported, kinematic viscosity and diffusivity set the relative thickness of thermal and kinematic boundary layers [measured by the Schmidt number $\text {Sc}$, defined in (2)], while the density difference produced by the scalar determines the driving force [measured by the Rayleigh number $\text {Ra}_{T}$, defined in (1)]. A key quantity to consider is the relative size of the pore space to the flow structure, which determines the penetration of the buoyant plumes, responsible of convective mixing, in the domain. The influence of these flow parameters on the convective transport efficiency, measured by $\text {Sh}$, is discussed here.

We initially consider a solid phase impermeable to the scalar, e.g. the case of solute convection. Accurate two-dimensional pore-scale simulations of Rayleigh–Bénard solutal convection are presented by Gasow et al. [31], where the porous medium is modelled as a matrix of aligned squares. They explored different values of porosity ($0.36\le \phi \le 0.56$) and Schmidt numbers ($\text {Sc}=1$ and $\text {Sc}=250$). The results, reported in Fig. 6 (green symbols) as measurements of Sherwood number, indicate fair agreement with two-dimensional Darcy simulations (black solid line, [54]). However, the pore-induced dispersion, which may be as strong as buoyancy, affects the flow structure and consequently $\text {Sh}$, and the scaling $\text {Sh}(\!\text {Ra})$ appears sublinear when the porosity is increased ($\phi =0.56$). At a low Schmidt numbers ($\text {Sc}= 1)$, pore-scale effects on the flow structure, e.g. wave number or width of the plumes, are qualitatively similar to those at high Schmidt numbers ($\text {Sc}= 250)$. In a complementary study, Gasow et al. [70] investigated large Schmidt numbers ($\text {Sc}= 250)$ convection, while focusing on the role of the medium properties. Results of [70] are reported in Fig. 6 (red symbols), and indicate that the dissolution coefficient depends on $\text {Ra}$ as

$$\begin{aligned} \text {Sh}=1+a\!\text {Ra}^{1-0.2\phi ^{2}} , \end{aligned}$$

(30)

where $a = 0.011 \pm 0.002$ is a pore-scale geometric parameter depending on shape and arrangement of the obstacles. The difference with respect to the Darcy case [simulations by [54]—black line, asymptotic best fit—Eq. (28)] is apparent, as it seems that within this range of parameters, systems with same $\text {Ra}$ (achieved with different values of porosity) exhibit very different convective transport properties.

Fig. 6

Sherwood number ($\text {Sh}$) as a function of Rayleigh–Darcy number ($\text {Ra}$) for two-dimensional pore-scale simulations. Results refer to solutal convection [31, 70] (i.e. solid impermeable to solute, green and red symbols) and thermal convection [30] (i.e. conductive medium, blue symbols). Results of two-dimensional Darcy simulations [54] (black solid line) and high-$\text {Ra}$ scaling [Eq. (28), grey solid line] are also shown

An additional degree of freedom is introduced by allowing a flux of scalar through the solid matrix, which may be the case for thermal convection. 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 control the penetration of the plumes into the boundary layer region, which in turn determines the heat or mass transfer rate. This physical mechanism has been described by three-dimensional pore-scale simulations of few pore spaces [71]. Later, in a complementary experimental study, Ataei-Dadavi et al. [72] 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 scaling derived for free fluids [73, 74] is recovered. The nature of this transition has been investigated in detail by Liu et al. [30] (blue symbols in Fig. 6). Within the frame of conductive media, two-dimensional direct numerical simulations have been used to investigate the microscale flow field at $\text {Sc}=4.3$. The obstacles consist of circles arranged in a regular manner. When the arrangement is not regular (not shown in Fig. 6), a slight decrease of $\text {Sh}$ is observed. In Fig. 6 it appears that the convective heat transport is less efficient compared to the configuration discussed before, in which the matrix was impermeable to solute (for a detailed discussion on the importance of the (im)permeability condition of the solid matrix, see [75]). In addition to the effect of the thermal conductivity of the solid, the measurements of [30] refer to relatively high values of porosity. As predicted by Eq. (30), the larger the porosity, the lower the $\text {Sh}$. The transition from porous convection to unconfined convection is controlled by two physical mechanisms, which are set by the properties of the porous matrix [30]. On the one hand, 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. On the other hand, the convection strength is reduced due the impedance of the obstacle array, leading to heat transfer reduction. The variation of $\text {Sh}$ with $\text {Ra}_{T}$ (not $\text {Ra}$) is reported in Fig. 7c, where the presence of these two distinct regimes is apparent. For sufficiently large $\text {Ra}_{T}$ or high porosity, the classical scaling is recovered ($\text {Ra}_{T}^{1/3}$, Grossmann and Lohse [73, 74]). When the Rayleigh number is lowered, however, the role of the porous structure in confining the flow is critical, and a correlation for the Sherwood number $\text {Sh}$ is proposed:

$$\begin{aligned} \text {Sh}\approx 1+c\phi \left( \frac{H}{\ell _{s}}\right) ^{4}\text {Sc}^{2}\text {Re}^{2}(\!\text {Ra}_{T})^{-1}, \end{aligned}$$

(31)

where the Reynolds number $\text {Re}$ is computed based on the velocity fluctuations and $c=8$ is a fitting parameter. This scaling is proved to be well approximated by $\text {Sh}\sim \!\text {Ra}_{T}^{0.65}$ [30]). The transition between these regimes appears clearly in Fig. 7d, when the compensated Sherwood number ($\text {Sh}/\!\text {Ra}_{T}^{-0.3}$) is shown as a function of the boundary layer thickness [$\delta =H/(2\text {Sh})$] divided by the average pore scale ($l_{s}$). The situation is schematically illustrated in Fig. 7a and b. When the thickness of the thermal boundary layer is comparable to the averaged pore length scale ($\delta /l_{s}=1$), the transition from one regime to the other occurs. 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 [76, 77].

Fig. 7

Pore-scale two-dimensional simulations [30]. a Exemplar dimensionless temperature field ($\vartheta$) (adapted with permission from [30]), being 0 and 1 the temperature values at the top and bottom boundaries, respectively. b Detail with explicit indication of the boundary layer thickness ($\delta =H/(2\text {Sh})$) and the average pore scale ($l_{s}$). The medium consist of aligned circular and conductive obstacles for Schmidt number $\text {Sc}=4.3$. c Sherwood number ($\text {Sh}$) is reported as a function of the Rayleigh number ($\text {Ra}_{T}$) for different values of porosity, $\phi$. Note that results for unconfined fluids ($\phi =1$) are also shown. d Compensated Sherwood number ($\text {Sh}/\!\text {Ra}_{T}^{-0.3}$) as a function of $\delta /l_{s}$

4 One-sided convection

The one-sided configuration introduced in Sect. 2.3 is representative of natural instances like geological CO$_2$ sequestration [12] and mixing in groundwater flows [78]. In these cases, a fluid-saturated porous domain (sketched in Fig. 4c-i) is allowed to exchange solute through the top boundary. The system is initially driven by diffusion [17, 60, 79]. The fluid layer below the upper boundary becomes progressively rich in solute, increasing the density of the liquid phase. When sufficiently thick, this high-density layer eventually becomes unstable and finger-like structures form [80, 81], evolve (i.e. grow and merge) and if the Rayleigh–Darcy number is sufficiently large [$\text {Ra}>O(10^3)$] the system may reach a quasi-steady regime. In this phase, the dimensionless solute flux $\widehat{F}$ computed as in Eq. (26), indicating the mass of solute dissolved through the top boundary per unit of surface area and time, is nearly constant over time. For simplicity, hereinafter we will refer to $\widehat{F}$ as the time-averaged value of flux in this constant flux phase. The role of the fingers in promoting solute mixing is crucial, as initially proposed by Ennis-King and Paterson [82] and Xu et al. [83], since the contribution of convection accelerates considerably the dissolution compared to the purely diffusive case. The domain progressively saturates with incoming solute, up to the point in which the local concentration difference between the upper fluid layer and the top boundary is reduced, and the dissolution rate suddenly drops. This phase is referred to as shutdown regime and it has been accurately described [14, 16, 17, 60, 84]. A thorough description of the whole dissolution process is provided by Slim [17].

In this section, we will review the results relative to Darcy and pore-scale flows in the one-sided configuration, and we will focus on the dependency of the dissolution rate $\widehat{F}$ on the flow parameters during the constant flux regime.

4.1 Darcy flows

When the Darcy model is considered [Eqs. (13)–(15)], the flow is uniquely controlled by the Rayleigh–Darcy number $\text {Ra}$, similarly to the Rayleigh–Bénard case discussed in Sect. 3.1. Numerical two-dimensional simulations agree on the value of the flux during the constant flux regime, which was initially determined by Hesse [56] to be

$$\begin{aligned} \widehat{F}=0.017. \end{aligned}$$

(32)

This observation has been later confirmed by a number of numerical studies [15,16,17, 60, 85,86,87]. We refer to [87] for a review of literature scaling laws in the presence of variations to this problem (anisotropy, geochemistry, etc.).

Fig. 8

Examples of one-sided studies. a Experiment with MEG in water in bead packs (adapted with permission from [19]). b Experiment with propylene glycol (PG) and water in Hele-Shaw cell (adapted with permission from [21]). c Darcy simulation with non-monotonic density profile (adapted with permission from [21])

In the instance of three-dimensional domains, the dynamics is analogue to that discussed above. However, due to the large computational costs, only few numerical works are available [15, 88, 89]. A seminal work in the field is presented by Pau et al. [15], who performed three-dimensional simulations and estimated the flux to be higher than in the corresponding two-dimensional case (32). These results refer to $\text {Ra}\le 9\times 10^{3}$, and additional data at larger Rayleigh–Darcy numbers are required to determine at the exact value for $\widehat{F}$, which has been estimated not to exceed 25% of the two-dimensional case [15, 87, 89]. It is apparent that the additional degree of freedom represented by the third spatial dimension adds significant complexity to the fingering phenomena [15, 88], with the flow structure being more complex and dynamical [58].

4.2 Pore-scale and Hele-Shaw flows

The determination of $\widehat{F}$ has been carried out beyond the Darcy model via pore-scale simulations and experiments, and via Hele-Shaw set-ups. As discussed in Sect. 3.1, a classical argument for Darcy convection requires that $\text {Sh}$ scales linearly with $\text {Ra}$ [67, 68]. The theoretical interpretation is that in natural convection $\text {Sh}$ is uniquely controlled by the diffusive boundary layer, and it is independent of the flow interior and any external length scale. Only for an exponent of one for $\text {Ra}$, i.e. $\text {Sh}\sim \!\text {Ra}$, it is possible to have an expression for $\text {Sh}$ that is independent of H [18, 68, 90]. As a result [see also Eq. (24)], the flux $\widehat{F}$ is expected to be independent of $\text {Ra}$, as it emerges from Darcy simulations [e.g. see correlation (32)]. Despite this robust theoretical framework, a different scaling for $\text {Sh}(\!\text {Ra})$ was found by many studies.

Most of the experimental studies investigating one-sided convection have been carried out with the aid of bead packs or Hele-Shaw cells, examples of which are reported in Fig. 8a and b, respectively. In the manner, the porous medium consists of a matrix of rigid spheres, typically made of a transparent material to allow optical access to the flow, enclosed in a transparent container. A Hele-Shaw cell, in turn, is obtained with two parallel and transparent plates separated by a narrow gap b (usually, less than 1 mm thick). When the fluid velocity in the cell is sufficiently low (gap-based Reynolds number $\ll 1$), the flow behaves as a laminar Poiseuille flow, i.e. the gap average velocity is proportional to the pressure gradient via the inverse of viscosity and to a constant equal to $k=b^{2}/12$, where k is defined as the equivalent permeability of the cell. Since this formulation represents an analogue of the Darcy law (14), the Hele-Shaw cell is commonly used as a tool to reproduce a flow through a porous medium. Bead packs and Hele-Shaw experiments used to derive scaling laws are discussed in the following.

A first $\text {Sh}(\!\text {Ra})$ scaling was proposed by Neufeld et al. [19], who used experiments in glass beads to mimic one-sided convection in porous media. The fluids employed were methanol and ethylene glycol (MEG) and water. MEG (upper fluid layer in Fig. 8a) is lighter than water when pure, but it presents a non-monotonic density profile as a function of the fraction of water. As a result, at the interface between the two fluids (identified by the white boundary in Fig. 8a), a heavier mixture forms and originates finger-like instabilities (white structures). The Sherwood number, estimated by tracking the receding interface between the two fluid layers, was measured to scale as $\text {Sh}\sim \!\text {Ra}^{4/5}$, and the result was explained with a phenomenological model based on a boundary layer theory: The lateral solute diffusion from the downward plumes into the upward ones is responsible for the reduction of local concentration gradients and the corresponding density differences driving the flow. This translates into a reduction of the flux, making $\text {Sh}$ to reduce with respect to the classical scaling. An analogue approach was employed by Backhaus et al. [20], who used Hele-Shaw cells and a layer of water located vertically above a layer of propylene glycol PG (a similar system is shown in Fig. 8b). They obtained the scaling $\text {Sh}\sim \!\text {Ra}^{0.76}$ and identified the plumes spacing as the key parameter controlling the Sherwood number. Similar results are derived by Tsai et al. [91] (Hele-Shaw and beads, $\text {Sh}\sim \!\text {Ra}^{0.84}$), Ecke and Backhaus [92] (Hele-Shaw, $\text {Sh}\sim \!\text {Ra}^{0.76}$) and Guo et al. [93] (Hele-Shaw, $\text {Sh}\sim \!\text {Ra}^{0.95}$). The discrepancy existing between these sublinear scalings and the linear theoretical [67, 68] and numerical findings in case of Darcy simulations [15,16,17, 60, 86] has been subject of active investigations.

To examine this mismatch, numerical simulations and theoretical arguments were used [21]. Accurate simulations were employed to mimic the behaviour of the fluids used in the experiments (characterised by a non-monotonic density–concentration curve, and with a concentration-dependent viscosity), which differ from the ones classically considered in Darcy simulations (linear dependency of density with concentration, constant viscosity). A snapshot of the concentration field obtained for a Darcy simulation with non-monotonic density profile is shown in Fig. 8c. It was found that the dissolution flux is determined by the mean scalar dissipation rate, $\widehat{\chi }$. Mixing in porous media has a universal character, and the nonlinear behaviour observed needs to be explained with effects not present in the classical Darcy–Boussinesq model. In particular, the authors observed that several differences exists between this simple Darcy model and the experiments reporting sublinear scalings. Among the others, they identified three main possible sources of discrepancies: (i) dependency of viscosity with the solute concentration, (ii) non-monotonic behaviour of fluid density with solute concentration, and (iii) compressibility effects (volume change during the process of dissolution). The conclusion of [21] is that while the concentration-dependent behaviour of viscosity has a minor effect, the role of the non-monotonic density–concentration profiles (shape of the density curves) may considerably affect the Sherwood number scaling law. The role of some of these fluid properties has been later investigated and will be discussed in the following.

The scaling analysis performed by Amooie et al. [90] for non-Boussinesq and compressible flows reveals that the scaling $\text {Sh}\approx 181.02+0.165\!\text {Ra}$ represents the best fitting for their data. Therefore, the authors propose that the previously reported sublinear relations could be in part a result of relatively limited parameter range of the simulations (as in the case of [94]) or in part because the Rayleigh–Darcy number of the experiments lies below the asymptotic limit, i.e. before the classical linear scaling establishes.

To avoid a non-monotonic dependency of density with concentration, i.e. to remove the fluid properties as a possible reason of nonlinear scaling, experiments in Hele-Shaw cells have been performed. Potassium permanganate (KMnO$_{4}$) and water are used as analogue fluids, with solid crystals of KMnO$_{4}$ placed on a metal grid located on top of the cell. Water gradually dissolves the crystals, which remain in a fixed position hold by the mesh, and the resulting interface between the light and the heavy fluid is always fixed and flat. This methodology, initially introduced by Slim et al. [79], allowed to cover a wide range of Rayleigh–Darcy numbers. In addition, variations of volume and fluid viscosity with solute concentration are negligible. Results by [95] report a linear scaling of $\text {Sh}$ with $\text {Ra}$. Later studies [25, 96] indicate that within the same value of permeability the scaling $\text {Sh}\sim \!\text {Ra}$ holds. In general, $\text {Sh}$ may still be a function of $\text {Ra}$ due to the presence of mechanical dispersion [97].

The works presented indicate that the fluid properties may not be sufficient to justify the nonlinear $\text {Sh}(\!\text {Ra})$ scaling observed. However, other physical mechanisms induced by the Hele-Shaw cell or the dispersion in the porous medium are not present in the classical Darcy model. These effects, labelled as finite-size effects, maybe be responsible of the nonlinear scaling observed, and will be discussed in detail in the Sect. 5.

5 Finite-size effects

Domain features like lateral confinement, thickness-induced Hele-Shaw dispersion and pore-scale dispersion have been identified to play a role on the nonlinear scaling of $\text {Sh}$ with $\text {Ra}$ or the flow structure. The influence of these finite-size effects on convection will be reviewed in this section.

5.1 Effect of confinement

A natural question arising from numerical simulations is what happens when the domain is confined in one of the wall-parallel directions, and we will address this topic here in the frame of Rayeligh–Bénard, the Rayleigh–Taylor, and the full reservoir-scale flow dynamics.

The flow in a porous Rayleigh–Bénard system at large $\text {Ra}$ consists of two distinct regions (see Sect. 3): (i) the near-wall region, characterised by the presence of protoplumes, and (ii) the interior of the flow, controlled by megaplumes. The average flow structure in each of these regions is quantified by via the time- and horizontally averaged wave number, k. While the near-wall region is hard to be described theoretically, the interior of the flow has been well characterised. In two dimensions, stability analysis [98] of the flow interior for $\text {Ra}\rightarrow \infty$ suggests that $k\sim \!\text {Ra}^{5/14}$, in fair agreement with numerical measurements that give $k\sim \!\text {Ra}^{0.4}$ [54]. In three dimensions, theoretical results [99] indicate that $k\sim \!\text {Ra}^{1/2}$, which is in excellent agreement with numerical measurements of [58, 64], who obtained $k\sim \!\text {Ra}^{0.52}$ and $k\sim \!\text {Ra}^{0.49}$, respectively. In addition, [59] observed with the aid of numerical simulations that supercells, representing clusters of protoplumes located near the boundaries, are the footprint of the megaplumes dominating the bulk of the flow. Unexpectedly, the correlation between these flow structures is observed to hold up to very high Rayleigh–Dacry numbers. This flow structure, however, may be considerably affected by the domain size.

Two-dimensional numerical simulations performed by Wen et al. [62] revealed that identifying the wave number may be complicated. Domains with low aspect ratio can dramatically reduce or even suppress convection. The study shows that the interior structure of a two-dimensional system may result strongly conditioned by the domain width, suggesting that the inter-plume spacing is not unique. The authors finally conclude that determining a precise high-$\text {Ra}$ scaling of the interior inter-plume spacing will require extremely long simulations in very wide computational domains.

Fig. 9

Influence of lateral confinement (domain width) on the development of the flow structures. Three-dimensional Rayleigh–Bénard simulations performed at $\text {Ra}=10^{4}$ are shown (adapted with permission from [58, 59]). a Dimensionless temperature distribution ($\vartheta$) in a cubic domain, being 0 and 1 the values at top and bottom boundaries, respectively, with gravity $\textbf{g}$ acting along z. Periodic boundary conditions are applied in the wall-parallel directions (x, y). The domain has dimensions L, W and H in directions x, y and z, respectively. The size of the domain is progressively increased in direction x, so that the aspect ratio increases from (panels b, c) to (panels h, i). For each value of , temperature fields taken at the centreline $(z=1/2)$ (c, e, g, i) and close to the bottom wall $(z=0.005)$ (b, d, f, h) are shown. Note that different colour bars apply to centreline and near-wall panels

In three dimensions, the effect of the domain confinement has been investigated by De Paoli et al. [58]. They performed numerical simulations at $\text {Ra}=10^{4}$ in Rayleigh–Bénard configuration, in domains having variable extension in one of the wall-parallel directions, namely x in Fig. 9a, and constant extension in the other directions ($W=H$). Periodic boundary conditions are applied in the wall-parallel directions. The relative size of the domain extension in directions y, z with respect to x is quantified by the aspect ratio . Four values of are considered in Fig. 9, with the domain progressively increasing in size from to . The corresponding temperature fields, taken at the centreline $(z=1/2)$ and close to the bottom wall $(z=0.005)$, are shown in Fig. 9b-(i). A strong confinement of the domain presents dramatic effects on the flow structures. For sufficiently large domains, e.g. , the near-wall cells reported in Fig. 9b are randomly oriented and show a wide distribution of sizes. When the domain width is progressively reduced, the cells are strongly constrained (Fig. 9f, h) and eventually end up in an extremely ordered pattern (Fig. 9d). The same applies to the flow structures at the centreline that for small domains () form sheet-like plumes. More quantitative results, estimated by means of the horizontal radial mean wave number of these simulations and additional larger domains (not shown here), indicate that the flow structures at the near-wall and in the interior of the flow are strongly constrained by the size of the domain. They found that at $\text {Ra}=10^4$ the flow is independent of the size of the domain for .

Decreasing the size of the computational domain in one direction will inevitably change the flow structure from a three-dimensional towards a two-dimensional character. This transition has been investigated in the frame of Rayleigh–Taylor instability by Borgnino et al. [100]. Among the other indicators, they analysed the evolution of the mixing length, i.e. the time-dependent vertical extension of the tip-to-rear finger distance, to determine whether the system exhibits a two- or three-dimensional behaviour. They observed that for sufficiently large Rayleigh–Darcy numbers ($\text {Ra}>10^{5}$), the growth of the mixing length is always linear in time in two and three dimensions (note that at lower Rayleigh–Darcy numbers the growth of the mixing length may be superlinear [57, 101]). The prefactor of the growth for the mixing length varies, being larger in two dimensions than in three dimensions. They performed three-dimensional numerical simulation with triply periodic boundary conditions, in which the dimension of the domain in a direction perpendicular to gravity, defined in the following “thickness”, is progressively reduced. Results indicate that when the thickness diminishes below a certain threshold value, the systems transitions from a three-dimensional to a two-dimensional behaviour. This critical value corresponds to the wavelength associated with the most unstable mode obtained from linear stability analysis [102, 103]. The sharp transition observed in this case is remarkably different than in turbulent convection [104]. In the turbulent case, the dimensional transition occurs dynamically, i.e. when the width of the mixing region exceeds the confined dimension, and it is smooth due to the coexistence of direct and inverse energy cascades.

The horizontal domain extension is also a parameter that dramatically affects the evolution of a buoyant current from injection to complete dissolution, e.g. in the configuration sketched in Fig. 1c relative to geological sequestration of carbon dioxide. Using the model for two-phase gravity currents proposed by MacMinn et al. [105], De Paoli [11] analysed the effect of the domain width on the maximum horizontal extension of the current of carbon dioxide. They performed two-dimensional simulations in which the domain width is progressively increased, while keeping the domain height and the volume of fluid injected constant. It was found that the layer of CO$_2$-rich solution may spread over a horizontal distance greater than 100 times the vertical extension of the layer, indicating that simulations are width-dependent, and very wide domains have to be considered (width to height ratio $\ge 140$).

5.2 Hele-Shaw flows

The working principle of the Hele-Shaw apparatus, briefly introduced in Sect. 4.2, is illustrated in Fig. 10a. The fluid is contained between two parallel plates separated by a narrow gap of thickness b, and the flow obtained in this configuration may be representative of a Darcy flow. When the flow is dominated by viscous forces (gap-based Reynolds number $\ll 1$), the depth-averaged fluid velocity is proportional to the vertical pressure gradient and to the inverse of the viscosity, in analogy to the Darcy law (14). This proportionality constant, defined as equivalent permeability of the cell, is $k=b^2/12$, and it used to draw a link between Darcy and Hele-Shaw flows.

In convective flows, the driving force of the system is the presence of a solute with concentration $C_0\le C\le C_1$, which produces a maximum density difference $\Delta \rho$ within the domain. In this frame, the analogy between Hele-Shaw and Darcy flow has been investigated quantitatively by Letelier et al. [97], who observed that a combination of fluid properties (Schmidt number, $\text {Sc}$), cell geometry (anisotropy ratio, $\epsilon =\sqrt{k}/H$) and flow velocity (U, defined in (8), which depends on $\text {Ra}$) determines the flow regime. They considered an incompressible flow (3), and averaged the Navier–Stokes and ADE equations, respectively (4) and (6), in the direction of the gap thickness to obtain the following dimensionless system

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

(33)

$$\begin{aligned}&\frac{\partial C^*}{\partial t^*}+\mathbf {u^*}\cdot \nabla C^*=\frac{1}{\!\text {Ra}}\nabla ^{2}C^*+\nonumber \\&\quad +\frac{2}{35}\epsilon ^{2}\!\text {Ra}\nabla \cdot \left[ (\mathbf {u^*}\cdot \nabla C^*)\mathbf {u^*}\right] , \end{aligned}$$

(34)

valid for $\epsilon$ small, $\text {Sc}\ge 1$ and $\epsilon ^{2}\!\text {Ra}\ll 1$. A linear dependency of density with concentration is considered. In this case $^*$ indicates dimensionless variables where the velocity scale is U defined as (8), the length scale is H, the timescale is H/U and the pressure scale is $\mu U H/k$. The concentration is made dimensionless as $C^* = (C-C_0)/(C_1-C_0)$ and $\textbf{k}$ is the unit vector with direction opposite to gravity. Equations (33)–(34) may be respectively interpreted as a Darcy law (14) and an advection–diffusion equation (15), both with additional corrective terms taking into account the contribution of inertia and solute redistribution due to the presence of the walls. In the frame of Hele-Shaw convection, three main regimes have been identified [97]: (i) Darcy regime (Fig. 10b) when $\epsilon \rightarrow 0$, the concentration profile across the cell gap is nearly uniform and the flow is well described by a Darcy model; (ii) Hele-Shaw regime (Fig. 10c) when $\epsilon \ll 1$, $\epsilon ^{2}\!\text {Ra}\ll 1$ and $\text {Sc}\ge 1$, characterised a gradient of concentration across the cell gap, but with one single finger; and (iii) three-dimensional regime (Fig. 10d), when the parameters do not fall in the above mentioned limits, the inertial effects become dominant and the fluid layer in the gap is unstable, so that multiple fingers appear across the cell thickness. It is apparent that the cell geometry plays a key role in determining the flow regime and that all laboratory experiments fall either in the Hele-Shaw regime or in the three-dimensional regime. With the aid of numerical simulations, Letelier et al. [97] provided an evidence for the reduction of the scaling exponent discussed in Sect. 4.2 for convective flows in the Hele-Shaw regime.

Fig. 10

a Front view of a convective flow in a Hele-Shaw cell in one-sided configuration [25], with the solute concentration being constant at top. Fluids consists of an aqueous solution of KMnO$_{4}$ (purple to black) and water (white). The reference frame (x, y, z) and the direction along which gravity ($\textbf{g}$) acts are also indicated. b–d Schematic representation of the side views of the cell (the thickness b is not to scale with respect to the height H). Three possible flow regimes as identified by Letelier et al. [97] are shown

These finding were later confirmed by the laboratory experiments of [25], where the flux has been measured for different values of permeability (i.e. different b). Note that when the Schmidt number is large (as in the case of [25], where $\text {Sc}= O(10^{3})$), the dispersive effects dominate over to the inertial terms. As a result, Eq. (33) reduces to the Darcy law (14) with additional dispersive corrections. These findings suggest within the Hele-Shaw regime the scaling exponent is affected by the anisotropy ratio $\epsilon$, as predicted by Letelier et al. [97], possibly explaining the discrepancy observed between Darcy simulations [21] and Hele-Shaw experiments [20].

Finally, the theoretical work proposed by Letelier et al. [97] has been recently generalised by Ulloa and Letelier [106] and Letelier et al. [107] to the case of more complex systems characterised by the presence of two layers of fluids with non-monotonic density profiles. The framework provided in [107] allows to evaluate and compare the mixing performance of different systems. They propose a universal law for the evolution of $\text {Sh}/\widehat{\chi }$, which is independent of the cell geometry ($\epsilon$) and directly proportional to $\text {Ra}$. Using this theoretical framework, they suggest that a possible reason for the sublinear scaling observed by Backhaus et al. [20] is the flow regime (Hele-Shaw regime) in which the experiments are performed.

5.3 Dispersion in bead packs

Recent developments in experimental techniques allowed accurate and non-invasive measurements of convective dissolution in three-dimensional porous media. The studies discussed in Sect. 4.2 are relative to thin domains, i.e. laboratory experiments in which the dimension of the cell in the direction perpendicular to the transparent walls is much smaller than the other two. This confinement may have an effect on the development of the flow structures (see Sect. 5.1) and on the dissolution efficiency of the system. We will present here three-dimensional measurements of convection in porous media, and discuss possible approaches to model dispersion in this context.

A remarkable contribution in the field on convection in three-dimensional porous media was presented by Lister [108]. This work is original because of the medium used, consisting of a fibrous material, and because of the remarkable visualisations performed. Beside this work, most of the investigations on convection in three-dimensional porous flows involved the presence of bead packs. The emergence of tomographic imaging systems over the last years has considerably sped up the research in this field. In a pioneering work by Shattuck et al. [109], magnetic resonance imaging (MRI) of three-dimensional convective flows in opaque media were presented, and plumes at low Rayleigh–Darcy numbers ($<20\pi ^{2}$) were visualised. Also X-ray computed tomography (CT) imaging scan is now frequently used to study mixing of miscible fluids. Wang et al. [110] and Nakanishi et al. [111] provided correlations for Sherwood as a function of Péclet and Rayleigh–Darcy number, and observed a sublinear scaling for $\text {Sh}$ with $\text {Ra}$, with exponent 0.40 and 0.93 respectively. The same methodology was employed by Liyanage et al. [112], who reported the emergence of characteristic patterns that closely resemble the dynamical flow structures produced by high-resolution numerical simulations. In a later study [113], the role of viscosity has been also investigated. While on the one hand [112, 113] observed that the flow is heavily influenced by dispersion, on the other hand a linear scaling $\text {Sh}\sim \!\text {Ra}$ holds, in contrast to previous studies. This discrepancy may be due to the relatively short range and small values of $\text {Ra}$ explored, which is well below the value in correspondence of which the system is observed to attain an asymptotic linear scaling [54, 59]. Employing the same measurement technique but different fluids, Eckel et al. [114] achieved larger Rayleigh–Darcy numbers ($\le$ 55,000). Through qualitative and quantitative observations of flow evolution, they also observed an enhanced longitudinal spreading of the solute, but in this case a sublinear scaling for $Sh(\!\text {Ra})$ holds.

These works agree upon the fact that dispersion is crucial in determining the $\text {Sh}(\!\text {Ra})$ scaling of the flow, and non-Darcy effects should be included in the models employed [115]. Dispersion has been identified as responsible for the early onset of convection [116]. In addition, Ghesmat et al. [87, 117] observed that the flow structures are influenced by the strength of dispersion and the dissolution rate $\widehat{F}$ is increased with increasing strength of dispersion. However, this finding does not apply in general and it seems to be limited to the range parameters considered [84]. With the aid of laboratory experiments, Liang et al. [26] proposed that, in addition to the Rayleigh–Darcy number, a flow with dispersion is controlled by a dispersive Rayleigh–Darcy number

$$\begin{aligned} \!\text {Ra}_{d}=\frac{UH}{\phi D_{\textrm{T}}}=\frac{\!\text {Ra}D}{D_{\textrm{T}}}, \end{aligned}$$

(35)

with $D_{\textrm{T}}$ the transverse dispersion, U the buoyancy velocity defined in Eq. (8) and H the domain height. In geological formations, assigning appropriate values to $D_{\textrm{T}}$ is not trivial, and it has been a debated topic (we refer to [48], for a thorough review on this subject). The anisotropy ratio $r=D_{\textrm{L}}/D_{\textrm{T}}$ (see Sect. 2.2.1) is also important to determine the flow character. As a result, the parameters space for convective porous media flows with dispersion is controlled by at least three parameters: $\text {Ra}$, $\text {Ra}_{d}$ and r. In order to quantify the relative importance of molecular diffusion to transverse dispersion, one can introduce the parameter [84, 118]

$$\begin{aligned} \Delta =\frac{ \!\text {Ra}_{d}}{\!\text {Ra}}=\frac{D}{D_{\textrm{T}}}, \end{aligned}$$

(36)

that can be used to rewrite the dispersion tensor (16) in dimensionless form as:

$$\begin{aligned} \frac{\textbf{D}}{D} = \textbf{I}+\frac{1}{\Delta U}\left[ (r-1)\frac{\textbf{u}\textbf{u}}{\vert \textbf{u}|}+\textbf{u}\textbf{I}\right] . \end{aligned}$$

(37)

This expression suggests that the case of pure diffusion is recovered when $D_{\textrm{T}}\ll D$, corresponding to $\Delta \gg 1$.

With specific reference to granular media, additional simplifications allow a further characterisation of the flow in the parameter’s space. Considering that the longitudinal dispersivity can be approximated [29, 118] as $\alpha _{\textrm{L}}=D_{\textrm{L}}/U\approx d$, we can rewrite Eq. (35) as

$$\begin{aligned} \!\text {Ra}_{d}=\frac{UH}{\phi D_{\textrm{T}}}=\frac{UH}{\phi U\alpha _{\textrm{T}}}=\frac{rH}{d}. \end{aligned}$$

(38)

In bead packs, the permeability can be inferred from the Kozeny–Carman correlation [45, 119], i.e.

$$\begin{aligned} k=\frac{d^2}{36k_C}\frac{\phi ^3}{(1-\phi )^2}, \end{aligned}$$

(39)

where $k_C=5$ is the Carman constant for monodispersed spheres randomly packed [120]. As a result, we can provide an expression for $\text {Ra}_{d}$ and $\text {Ra}$ that is an explicit function of the domain (g, H), fluid ($D,\mu ,\Delta \rho$) and medium ($\phi ,d,r$) properties. This information is particularly important when we characterise the flow in the three-dimensional parameters space $(\!\text {Ra}_{d},\!\text {Ra},r)$, which we will do in the following.

Fig. 11

Parameters space $(\!\text {Ra},\!\text {Ra}_{d})$ with indication of iso-$\Delta$ lines (blue lines). With respect to experiments in bead packs, if only d varies, the parameters $(\!\text {Ra},\!\text {Ra}_{d})$ are locked to the green curves. An example for a realistic range of parameters is shown by symbols (circles), where each line corresponds to one value of density contrast ($\Delta \rho$). Conversely, if only $\Delta \rho$ is varied in the experiments, $(\!\text {Ra},\!\text {Ra}_{d})$ are locked to horizontal lines (diamonds). Triangles indicate the slope of the green and blue lines

First, we reduce the parameters space to $(\!\text {Ra},\!\text {Ra}_{d})$ by taking into account [48] that $r=O(10)$ may be a reasonable approximation for advection-dominated systems, and we consider $r=10$. Note that no remarkable difference in the flow structure and $\text {Sh}$ occurs for $r>10$, provided that $\text {Ra}$ and $\text {Ra}_{d}$ are sufficiently large (namely, $10^{4}$ and $10^{3}$, respectively [84]). For $r\le 1$, the flow is qualitatively similar or that observed in the absence of dispersion [54]. With respect to the remaining parameters, $\text {Ra}$ and $\text {Ra}_{d}$, we can rewrite both as a function of the beads diameter and find that $\text {Ra}_{d}\sim 1/d$ and $\text {Ra}\sim d^{2}$. This implies that if we consider an experiment in which only d varies, we are locked to one of the green lines of the parameters space $(\!\text {Ra},\!\text {Ra}_{d})$ shown in Fig. 11, corresponding to $\text {Ra}_{d}\sim \!\text {Ra}^{-1/2}$. Using realistic laboratory properties, we obtain that a possible range for the experimental parameters $(\!\text {Ra}_{d},\!\text {Ra})$ at variable d consists of the circles in Fig. 11 (each series of circles corresponds to one value of density difference, $\Delta \rho$). Alternatively, we consider the case in which the medium is fixed (d constant) and the fluid density contrast varies. Since $\text {Ra}_{d}$ is independent of any fluid property, a variation of $\Delta \rho$ will correspond to a horizontal line of the parameters space (red symbols in Fig. 11, in which each series is a different d). Finally, we consider the case of a constant value of $\Delta$. It follows that this is achieved when $(\Delta \rho )^{-1}d^{-3}$ is constant [blue lines in Fig. 11]. This condition is extremely challenging to be obtained experimentally because it implies a simultaneous variation of $\Delta \rho$ and d.

Fig. 12

Concentration distribution at high $\text {Ra}$ and $r=10$ (Rayleigh–Darcy number and $\Delta$ indicated within each panel). a and c Rayleigh–Bénard configuration (adapted with permission from [84]). b and d One-sided configuration (adapted with permission from [26]). When $\Delta \gg 1$ (a, b), plumes grow vertically in a symmetric fashion (columnar flow). When $\Delta \ll 1$ (c, d), dispersion makes the plumes to expand in the horizontal direction (fan flow)

With the aid of numerical simulations the problem of decoupling two of the governing flow parameters $(\!\text {Ra},\!\text {Ra}_{d})$ can be solved, and their relative effect on the $\text {Sh}$ or $\widehat{F}$ can be investigated. Wen et al. [84] considered a Rayleigh–Bénard configuration and investigated systematically a range of flow parameters indicated in Fig. 13 (red squares). The flow structure is mainly ruled by $\Delta$, which determines the mechanism controlling convection. If $\Delta >10^{5}$, the flow is ruled by molecular diffusion, plumes grow symmetrically (Fig. 12a) and the structure is analogue to the symmetric flow observed in the absence of dispersion (see Fig. 4b-i). When $\Delta <1$, mechanical dispersion dominates over convection, and its inherent anisotropy ($r\gg 1$) sets the non-symmetric flow structure (fan flow) shown in Fig. 12c, in which plumes widen as they move away from the wall. A similar behaviour is observed in the corresponding one-sided cases (Fig. 12b, d) by Liang et al. [26].

The effect of the flow structure on the Sherwood number was also quantified. Note that in case of dispersive flows the Sherwood number contains the magnitude of the velocity at the top wall in its definition. Indeed, while the vertical component of velocity in zero at top (no penetration), a nonzero velocity parallel to the wall is admitted (free-slip), which produces solute spread due to dispersion. Alternatively, $\text {Sh}$ can be inferred from the time derivative of the total mass of solute in the domain. We report in Fig. 13 a visual interpretation of the dominant mechanism in each region of the $(\!\text {Ra},\!\text {Ra}_{d})$ space, where regions controlled by different mechanisms are separated by dashed lines. Liang et al. [26] found that in Rayleigh–Bénard configuration, when $\Delta >O(1)$ molecular diffusion dominates over mechanical dispersion, although a small contribution of mechanical dispersion may increase $\text {Sh}$. When $0.02<\Delta <O(1)$, both mechanical dispersion and molecular diffusion determine the value of $\text {Sh}$. A linear scaling $\text {Sh}\sim \!\text {Ra}$ holds when $\Delta <0.02$, but $\text {Ra}_{d}$ is also important since it determines the prefactor of the scaling law, as it has been later observed by Erfani et al. [121].

These results for dispersion-dominated flows ($\Delta <0.02$) could also provide an additional interpretation to the Hele-Shaw experiments of [25], where within one value of cell gap (i.e. one value of permeability and mechanical dispersion), the flux remains nearly constant, i.e. the prefactor is constant. On the other hand, Hele-Shaw flows do not exhibit transverse dispersion [122], therefore we believe that such one-to-one comparison between results of dispersion for Hele-Shaw and bead pack flows may not be appropriate.

Fig. 13

Range of parameters space explored with simulations [84, 117, 118] and experiments [19, 26] with dispersion. The configuration (one-sided—OS or Rayleigh–Bénard—RB) is indicated. All cases refer to $r=10$, with the exception of [117] where also different values of r are considered. Effect of r on convection is also discussed in [84], but the corresponding data are not reported in this figure. In this parameters space, $\Delta$ sets the flow behaviour: diffusion-dominated [$\Delta <O(1)$], dispersion-dominated [$\Delta >0.02$] or influenced by both diffusion and dispersion [$0.02<\Delta <O(1)$]

Recent works investigated the role of mechanical dispersion with the aid of simulations. A possible complementary approach with respect to the formulation used by Wen et al. [84] and Liang et al. [26] is proposed by Dhar et al. [123], where the dispersive Rayleigh number is replaced by a parameter quantifying the strength of longitudinal dispersion compared to molecular diffusion. More recently, Tsinober et al. [118] performed simulations in one-sided configuration. They modelled fluids with constant viscosity and linear density–concentration profiles, and derived a linear correlation between $\text {Sh}$ and $\text {Ra}$, where the prefactor is a function of molecular to dispersive Rayleigh–Darcy numbers, $\text {Ra}/\!\text {Ra}_{d}$. This correlation fits well their results, but it does not fully capture the trend predicted by Wen et al. [84]. This difference may possibly be due to several reasons, including the parameters space (and perhaps the different regimes) explored compared to [84] as it appears from Fig. 13, where the parameters investigated in some of these studies are reported. Each of these works involves a specific configuration (Rayleigh–Bénard or one-sided), a different formulation (different model for the fluids and different dimensionless parameters) and a different region of the parameter space. Therefore, providing a precise and general description of convective and dispersive flows in porous media is still not possible, and further studies systematically investigating a broad range of the $(\!\text {Ra},\!\text {Ra}_{d})$ space are required.

Finally, an novel approach consists of including the effects of dispersion also in the momentum equation. Gasow et al. [124] used two-dimensional pore-scale and Darcy simulations to study a Rayleigh–Bénard flow. 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, proportional to $\sqrt{k}$, and the domain size, H). The authors proposed a two-length-scale diffusion model, in which the pore-scale dispersion is accounted into the momentum transport as a macroscopic diffusion term. A similar model, which is found to be valid for a wide range of porosity values and is based on the effective viscosity, has been proposed to account for pore-scale effects in advection-dominated systems in the absence of convection [125].

6 Summary and future perspectives

In this work, we have reviewed recent developments on convection in porous media. We focused on state-of-the-art measurements of dissolution and mixing in archetypal flow configurations. Despite the well-known mathematical formulation of the problem, the role that several physical processes (e.g. finite-size effects) have on the dissolution and mixing is not yet fully understood. This is also due to the great complexity of the physics involved: convection in porous media is a nonlinear phenomenon taking place in multiphase and multiscale systems, eventually located thousands of metres beneath the Earth surface. Notwithstanding the intrinsic difficulties associated with performing reliable measurements in such systems, remarkable developments have been achieved in recent years.

The porous Rayleigh–Bénard configurations, consisting of a fluid-saturated porous slab with fixed density at top and bottom boundary, has been extensively investigated [53, 54, 62, 63]. The governing parameter of the flow is the Rayleigh–Darcy number $\text {Ra}$, a measure of strength of convection relative to diffusion. Three-dimensional Darcy simulations performed at unprecedented Rayleigh–Darcy numbers, $O(10^{5})$ have been used, and the existence of new flow features labelled as supercells emerged [58, 59]. Two-dimensional and three-dimensional simulations have shown that ultimately a linear scaling of the dimensionless dissolution coefficient is attained, namely $\text {Sh}\sim \!\text {Ra}$. While in the two-dimensional case [54, 62] this scaling sets in at $\text {Ra}\le 10^{4}$, in three-dimensional flows [58, 59] the ultimate state is expected to take place at $\text {Ra}\ge 5\times 10^{5}$, which is beyond the present numerical capabilities. Pore-scale simulations have revealed a more complex scenario, in which the heat/mass transfer is also influenced by porosity [30, 31], Schmidt number and relative conductivity of fluid and solid phases [76, 77]. These extensive numerical campaigns have led to the development of physics-based correlations for $\text {Sh}$ as a function of the flow parameters. In addition, the relative size of boundary layer and average pore space has been identified as a critical flow feature controlling pore-scale convection [30].

The second archetypal configuration considered is the one-sided configuration [16], where solute dissolves in an initially solute-free porous domain from the upper boundary, with all other boundaries being impermeable to fluid and solute. The flow is characterised by an intermediate phase in which the dissolution rate $\widehat{F}$ is quasi-steady. While Darcy simulations report a constant $\text {Ra}$-independent value $\widehat{F}$ [16, 17, 56, 60], experiments in bead packs [19] and Hele-Shaw cells [20] revealed that $\widehat{F}$ is a function of $\text {Ra}$. The discrepancy observed has been attributed to non-Darcy effects present in the experiments and not accounted by the simulations [21]. This has stimulated further studies focusing on the role of finite-size effects observed in Hele-Shaw [25, 97, 106, 107] and bead packs experiments [26, 84]. The analysis of recent numerical and experimental results [26, 84, 118] highlights the complexity of this system, which is controlled by at least three parameters, respectively quantifying the relative strength of (i) convection and diffusion ($\text {Ra}$), (ii) convection and dispersion ($\text {Ra}_{d}$), and (iii) longitudinal and transverse dispersion (r). The huge parameters space defined in this way and the need for both numerical and computational studies represents a major challenge in this field.

Improvement of numerical and experimental techniques allowed a detailed characterisation of the flow and a better understating of the phenomena involved. The combination of theoretical modelling, numerical simulations and laboratory observations will pave the way to derive and validate large-scale models to be employed in real geophysical and engineering situations. These findings will be crucial to tackle problems associated with grand societal challenges, such as energy transition and climate change mitigation [7].

To conclude, in Sect. 6.1 we will briefly review recent advancements in experimental techniques, and in Sect. 6.2 we will also discuss the importance of additional effects not considered in previous sections of this paper.

6.1 Recent developments in experimental techniques

One intrinsic challenge associated with measurements in porous media consists of the impossibility of optically accessing inner regions of the flow. An overview of the experimental techniques available to perform measurements in opaque media is presented by Poelma [23]. Among the different imaging techniques employed for porous media [36], magnetic resonance imaging (MRI) [126,127,128] and X-ray tomography [110, 129] are the most common, and allow to obtain non-invasive and non-intrusive three-dimensional measurements of inner flow regions. Despite the advantages mentioned, these techniques are high-priced and typically lack in resolution in both space and time, making fast and small-scale flows hard to measure. However, thanks to the recent technological progresses, these measurement techniques allowed a detailed characterisation of both medium and flow also at small scales. Some examples are the X-ray synchrotron microtomography [130], with resolution in space of 3.25 $\upmu$m and in time of 6 s. Recent experiments [131, 132] have shown that the resolution can be further lowered down to 2.3 $\upmu$m, with a technique also allowing for higher resolution in time. At the time being, similar performances are also achieved by commercial micro-CT systems. Optical measurement in three-dimensional porous media can be also performed by matching the refractive index of fluid and medium [115, 133, 134], provided that a suitable fluid is available. This is not always granted, since fluids with refractive indexes of interest may come with side effects such as high costs or high hazard [135].

Additional challenges associated with laboratory experiments, in particular with respect to geological sequestration of carbon dioxide, consist of reproducing realistic porous media and ambient conditions. For instance, at the depths at which CO$_{2}$ is supercritical, the pressure is of the order of tens of bars, and performing controlled experiments with optical access is not trivial. This obstacle has been recently successfully overcome [61, 136], and the methods proposed may represent a first important step towards investigations in more complex geometries. With respect to the design and production of synthetic media at the laboratory scale, microfluidic devices mimicking porous materials are usually made of polydimethylsiloxane (PDMS), which has the drawback of being permeable to CO$_{2}$. A solution has been recently proposed by De et al. [137], who developed a new method to fabricate a two-dimensional porous medium (regular array of cylinders), consisting of bonding of a patterned photo-lithographed layer on a flat base. Additional examples of manufacturing techniques for analogue porous media are provided in [138]. Real geological formations are inherently disordered and heterogeneous, and mimicking this feature in laboratory models is essential to capture the role of the medium heterogeneities on the solute mixing. The technique proposed by Guo et al. [139] addresses this issue, and it consists of a cell made of 3D printed elementary blocks designed to be easily rearranged to obtain a desired permeability field.

Finally, we conclude with an overview of recent developments in experimental techniques employed in Hele-Shaw cells. The relative low cost and ease of implementation makes this apparatus widely employed to study buoyancy-driven flows. Classical optical methods based on light intensity measurements of patterns induced by density (or density gradients) fields, such as Schlieren and related techniques [24, 140], have been combined or improved to increase the accuracy of the measurements performed. Accurate temperature [141, 142] and concentration [79, 95, 96] measurement techniques have been recently introduced. Velocity measurements have been also performed using advanced particle image velocimetry (PIV) and particle tracking velocimetry (PTV) techniques specifically designed for Hele-Shaw flows [143, 144], or with the aid of machine learning techniques, namely convolutional neural network (CNN) [145]. A separate (i.e. not simultaneous) measurement of scalar and velocity fields complicates the analysis of the phenomena involved and the description the underlying physical mechanisms. Recently, novel techniques for simultaneous temperature/concentration/velocity measurements have been proposed [96, 146], which are particularly useful to enable reliable comparisons against numerical findings.

6.2 Additional effects influencing mixing

Convection and mixing in real engineering and geophysical problems are far more complex than the idealised conditions depicted in this review, due to the non-ideal medium, fluid, and ambient conditions. Here we will discuss the influence of conditions not present in the configurations previously discussed, and we will provide some references for the interested readers.

We focused on processes in which the Boussinesq approximation applies, i.e. the density variations induced by the presence of a scalar are only significant within the gravitational term of the momentum (Darcy) equation, and can be neglected elsewhere. In general, this may not be the case, and a criterion for the applicability of the Boussinesq approximation has been derived [28]. For instance, in case of iso-thermal brine transport, fluid volume changes may be neglected when $\text {Ra}\tilde{\rho }_{r}/\Delta \rho \gg 1$, being $\Delta \rho$ the maximum density difference and $\tilde{\rho }_{r}$ the reference fluid density. Interestingly, this condition is independent of $\Delta \rho$ and it is widely fulfilled for geothermal processes, when $\text {Ra}\approx 10^{1}-10^{3}$ and ${\rho }_{r}/\Delta \rho \approx 10^{2}-10^{3}$ [147]. Numerical simulations of the fully compressible CO$_{2}$ sequestration process suggest that compressibility and non-Boussinesq effects do not significantly impact spreading and mixing [90]. An aspect particularly relevant when considering experiments with analogue fluids is that the mixing rate strongly depends on the shape of the fluids density–concentration curve and, in particular, on the position of the maximum of this curve [21]. This effect, along with volume variations in the fluid phase [148], may influence the dynamics of the mixing process, and the findings discussed in this review cannot be generalised to fluids with a non-monotonic density–concentration profile or in the presence of significant volume variations.

Geological formations are typically characterised by anisotropic and heterogeneous media. The effect of anisotropy has been well characterised by assuming that the permeability tensor is anisotropic [149, 150], and it has been shown that anisotropy is in general favourable since it increases the rate of dissolution and anticipates the onset of convection [60, 63, 151]. These studies assume that a preferential direction exists, i.e. the permeability tensor takes a diagonal form in a reference frame that usually has a direction aligned with gravity. This simplified model does not take into account that formations have heterogeneities, which are also source of anisotropy, and discerning these two features of the medium represents a strong simplification. It has been proposed [89, 152] that the model of anisotropic medium discussed above (in which the permeability tensor is diagonal in some reference frame) may represent a good candidate to investigate heterogeneous media. Different models for heterogeneous formations have been introduced, consisting of essentially three categories: spatially variable permeability fields [153] (with no preferential direction), long and thin impermeable barriers [89, 152, 154], and layered formations (i.e. regions in which high- and low-permeability strata alternate) [155,156,157]. Although a general model for convection in heterogeneous media is not available yet, these studies provide an initial framework to understand the long-term behaviour of these systems.

The role of the fluid properties may also affect the flow evolution and solute mixing. The effect of viscosity, for instance, may be crucial in determining the stability of a layer, and we refer to [158, 159] for a review on this topic. Another effect that is increasingly studied is the reactivity of the medium with the fluid: the solute present in the fluid may induce dissolution or precipitation, which corresponds to a variation of the medium porosity and permeability. Recently this problem has been actively investigated [121, 160, 161], also due to the improvement of numerical capabilities. It has been reported [87] that medium morphology modifications occurring in the presence of convective flows affect solute mixing in non-trivial manners. Cardoso and Andres [162] showed that the reacting system rock-CO$_{2}$ may be described by a first-order chemical reaction stimulating numerous studies on convective-reactive porous media flows, reviewed in [163, 164].

Finally, the effect of the ambient flow conditions may be also important [165]. It was observed that the presence of a background flow influences the onset of convection [166]. Experiments in one-sided configuration [167] revealed that while convection may be hindered and suppressed, dispersion enhances, with an overall contribution with respect to flux in the absence of background flow that can be positive, negative or neutral. Tsinober et al. [168] observed with the aid of simulations that three regimes exist, in which convection dominates, background flow dominates, or these two contributions have the same strength. These results are relevant, since they can contribute to derive new models suitable for prediction of dissolution at the scale of the reservoir [11, 105, 169] and through the entire lifetime of a buoyant current in a porous formation [169,170,171,172].