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


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 densitydriven 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. 1(a), approximately a meter 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 buoyancyinduced 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].
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. 1(c), 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. 1(c), 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.
trategy involved adding a third loride-to a solution of water and amn that case, the added component inthe fluid released upon solidification.mptly shut off the convective instaformation of unwanted freckles.an extra component to a metallurgie forming a solid product easier, the sideration is whether that product properties.Thus, to make quantitaulticomponent alloys of metallurgit address their increasingly complex at need led metallurgists at NIST to with researchers at the University of freckle-predictor criterion, a Rayleigh ure for predicting whether freckles t nickel-based superalloys. 13Because ulticomponent systems are typically uum models 14 based on conservation , energy, and species, all coupled to ns.The models are amenable to come the one shown in figure 3.

s and distinct mushy layers
s ternary systems constitutes a small but rd understanding complex multicompot studies of systems with relatively simple e revealed new, and in some cases unexssociated with convection in mushy layers.stems, the first ingredient for understandsolidification is the ternary phase diagram.4a is for a system with components A, B, A represents a solvent, such as water, and , such as salts.Upon cooling from below, A solidifies at point P to form single-phase dendrites in a primary mushy layer, shown in figure 4b, with B and C le in the liquid.The liquid mass fractions of both B and C increase with depth in the primary mushy layer as A is removed from the liquid phase to form the dendrites.At point S, components A and B solidify together, forming two-phase dendrites in a secondary mushy layer, with component C rejected.FIGURE 3. SEA ICE GROWTH is simulated here using an enthalpy method in an effectively two-dimensional Hele-Shaw cell of depth H.The cell is cooled from above. 15Salt-enriched plumes (yellow) with dimensionless salinity S drain from brine channels in sea ice with porosity χ into the underlying liquid (blue).

R CONVECTION
ushy layer with cold, elow warm, salty water, t diagram.The layer is d solutally unstable, a supereutectic growth.ent is water, or solvent pon solidification; the is salt, or solute B, and ssume that heat difhan solute B. isplaced fluid parcel sponse to its new enntain thermodynamic ieve the composition uidus constraint, the e of the surrounding lution may lead to a e accompanying per-, although secondary n erosion, generates of the flow into nars.Now consider a ternary mixture whose primary mushy layer is thermally and solutally stable-specifically, the bottom is cold and rich in two salts, solutes B and C, as shown in the right diagram.Assume heat diffuses much faster than solutes and that solute C diffuses much faster than solute B. For simplicity, suppose that buoyancy forces act only with respect to solute C.
An upwardly displaced fluid parcel warms rapidly but remains relatively rich in the slowest-diffusing component, solute B. Thermodynamic equilibrium-the ternary liquidus constraint-is maintained by a relative depletion of solute C through solutal diffusion.The parcel is thus fresher in solute C and lighter than its local environment, so it experiences an upward buoyancy force that drives an instability, despite the apparently stable background stratification.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 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 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 setups are not captured by the classical Darcy formulation [21].

Cold, fresh
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 Sec. 2 we describe the mathematical models and the idealized configurations used to investigate convection in porous media, and we derive a unified formulation to evaluate and relate mixing in different flow configurations.In Secs. 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 Sec. 5. Finally, in Sec.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 Secs. 3 and 4.

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, ∆ρ, 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 Ra T defined on the characteristic length scale of the flow H where g is the acceleration due to gravity, D is the molecular diffusivity and µ is the fluid dynamic viscosity.The ratio of kinematic viscosity to solute diffusivity D (or molecular diffusivity) determines the Schmidt number with ρ r the average (or reference) fluid density within the domain.Similarly, for thermally-driven flows one can define the Prandtl number (Pr ), in which the molecular diffusion is replaced by its thermal counterpart.Modelling heat or mass transport at the porescale requires to resolve the flow within the interstitial space.Momentum transport is controlled by continuity and Navier-Stokes equations, respectively: where ũ, ρ and p are the velocity, density and pressure fields, respectively, and 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 Sec.6.2).The fluid density ρ 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 linearized, the EOS may be rewritten to obtain the density ρ as a function of temperature ( T ) and concentration ( C) (possible limitations of this approach are discussed in Sec.6.2).With respect to the value ρr defined at the reference state ( Cr , Tr ), it reads where α c , α t are the expansion coefficients relating the density to the variations of concentration and temperature, respectively, being typically α c > 0 and α t < 0. In this case, assuming the presence of solute scalars only, Eq. ( 5) reduces to the form ρ = ρ( C) = ρr + α c ( C − Cr ).Solute conservation is accounted by the advection-diffusion equation: Eq. ( 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 Sec.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.

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 [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.
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 ϕ, 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: 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 ∆ρ, which is defined as We observe that U is independent of the any length scale, and it relates to the fluid (∆ρ, µ), medium (k) and domain (g) properties.In addition to the domain length scale (H), one can consider as a reference length scale the distance ℓ over which advection and diffusion balance [17] (in thermal convection, ℓ = 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 Ra obtained combining the Rayleigh number Ra T (1) and the Darcy number Da (7).In particular, Ra = Ra T Da /ϕ in the instance of solutal convection, with the solid being impermeable to the solute fluxes, and Ra = Ra T Da for thermallydriven 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 solutefree.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, ϕ) and fluid (D, µ, ∆ρ) properties.In addition, when the spatial coordinates are made dimensionless with respect to ℓ (9), 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 √ k.In quantitative terms, the criterion above is fulfilled when: (i) the pore-scale Reynolds number is small, i.e., viscous dissipation (µU ) dominates over inertia (ρ r U 2 √ k), and (ii) the smallest length scale of the flow (ℓ) is large compared to the pore size ( √ k).These conditions are matched if: i.e. when Reynolds (Re) and Péclet (Pe) numbers are much less than unity.Note that in these definitions the pore length scale ( √ 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 ϕ and permeability k (Fig. 2a) fulfilling the conditions (11)- (12).The flow field is fully described by the continuity and Darcy equations, respectively: Note that in this case 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 (ũ) via the Dupuit-Forchheimer relationship u = ϕũ [18].Same applies to pressure p and density ρ.
The evolution of the concentration field is controlled by the advection-diffusion equation 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 Sec.2.2.1 this coefficient may be replaced by a dispersion tensor D, that depends on the local flow conditions (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.

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 homogenizing 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 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 L , longitudinal dispersion) and perpendicular (D T , transverse dispersion) to the local flow velocity u.This model is labelled as Fickian dispersion [38].With these assumptions, the dispersion tensor takes the form: where I is the identity tensor, and the coefficients α L = D L /U and α T = D T /U correspond to the dispersivities of the medium in longitudinal and transverse directions, respectively.For solute transport and Pe ≫ 1, dispersion in the cross-flow direction is typically 1 order of magnitude smaller than in the stream flow direction (possible limitation of this assumption are discussed in [39]).The ratio of these two contributions is quantified by the dispersivity ratio The magnitude of D L and D 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], 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 setups 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 Sc ≤ 550.For a The correlations proposed by [42] have been employed.The advective flow is divided in several regimes, discussed in the text.
monodispersed close random packing the porosity is ϕ = 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 τ = 0.68 [see 47, and references therein].We use the empirical correlations proposed by [42], obtained for laboratory experiments, i.e., in architecturecontrolled media, whereas we refer to [48] for a review of dispersion relations in field-scale data.The results proposed by [42] are valid for liquids and at Sc ≤ 550, and we report the dispersivity ratio in Fig. 3 (for Sc > 550, similar correlations are provided).Four main flow regimes have been identified, for increasing P e: (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; (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 therin]) Péclet numbers.A similar regime classification has been also proposed by [50].
With this example we have shown that in general r varies with P e 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 L and D T or r [26], which is a reasonable approximation if a narrow range of P e is considered.Results on the effect of dispersion on convective flow are presented in Sec.5.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 analyze the evolution of the system.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. 4(a).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 (u • n = 0, being u the fluid velocity and 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 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 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: where ⟨•⟩ indicates the volume average.Eq. ( 18) relates the mean squared concentration, the solute flux through the walls F and the mean scalar dissipation within the domain χ, respectively defined as with L domain width, and 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.Eq. ( 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 (χ is always positive, therefore it contributes to a reduction of scalar variance ⟨C 2 ⟩).
To enable comparisons among different systems, a possible set of dimensionless variables consists of L for lengths, ϕL/U for time, and U for velocities.The concentration C is made dimensionless as C = (C − C 0 )/∆C, where • indicates dimensionless quantities and ∆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: with F = F 1 L/(D∆C) the dimensionless flux and χ = χL 2 /[D(∆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 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 where • indicates the time-averaging operator.We observe in Fig. 4(b-ii) that while a non-zero instantaneous contribution d⟨ C 2 ⟩/d t is present, F and χ fluctuate around their timeaveraged value (black dashed line).Note that the reference length scale normally used in this configuration is L = H, which gives in (21) the prefactor ϕD/(U L) = 1/Ra.The quantity used to evaluate the mass transfer in this configurations is the Sherwood number defined as the relative contribution of convective and diffusive to diffusive mass transport.
Using the definition of F and Eq. ( 22), Sh can be related to the flux and the mean scalar dissipation [53]: ii. One-sided (Fig. 4c-ii): the domain is impermeable to solute at the lower wall (F 0 = 0).By setting L = ℓ as defined in (9), Eq. ( 21) is independent of Ra and reads where This choice for L is convenient to compare systems having different Ra because the value of F appears to be universal, as will be later discussed in Sec. 4. The time-dependent flow originated from this configuration consists of three-main flow regimes [17,60,61].Initially ( t < 10 3 ) diffusion dominates and a highconcentration high-density unstable fluid layer thickens.At a later stage ( t < 16Ra), convection takes place and plumes formed at the top boundary layer grow and invade the domain.
In this phase F is statistically steady and characterized by a value (black dashed line) that is independent of the Rayleigh-Darcy number considered.A similar behaviour holds for χ, but a closer inspection reveals that after the fingers reach the bottom ( t > 10Ra) an increase of d⟨ C 2 ⟩/d t is observed.A corresponding decreasing behaviour is reflected in χ, but with half the amplitude.After the upper layer of the domain is also saturated ( t > 16Ra), the dissolution flux F drops, and the system enters the shutdown regime.iii.Rayleigh-Taylor (Fig. 4d-ii): the domain is impermeable to the solute ( F = 0) and Eq. ( 21) reads 1 2 The flow is initialised considering two fluid layers of different density in an unstable configuration.Eq. ( 27) suggests that all the potential energy initially stored by keeping the two phases segregated is dissipated as time evolves.Both ℓ and H can be considered as reference length scales, depending on which part of the flow evolution is considered.However, [52] have shown that L = ℓ provides a universal picture for the evolution of χ, and results are presented in Fig. 4(d-ii) using this length-scale.Similarly to what observed in the one-sided configuration, the flow is initially controlled by diffusion ( t < 10 3 ).Afterwards (10 3 < t < Ra /2) the formation of fingers is observed, which merge and grow, accelerating mixing.In this phase, occurring at Ra < t < 3Ra in the simulations considered, χ is observed to increase in Darcy simulations, as shown in Fig. 4(d-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, 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 ( F , Sh or χ) with the flow parameter (Ra).These results will be reviewed Secs. 3 and 4 for the Rayleigh-Bénard and the one-sided configurations, respectively.

Rayleigh-Bénard convection
Rayleigh-Bénard convection produces the statistically-steady flow discussed in Sec.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.

Darcy flow
In the Darcy case [Eqs.( 13)-( 15)], the system is uniquely controlled by the Rayleigh-Darcy number Ra, defined in Eq. (10), which sets the flow structure.The behaviour of Sh withRa 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 Ra < 4π 2 , the mass transport is purely diffusive [65,66] and no convective motion arises (Sh = 1).The flow is maintained quiescent by the dissipative (diffusive) effects that dominate over convection.For increasing Ra, instabilities appear in the form of steady convective rolls [55] with corresponding increase of the convective mass transfer.When Ra ≈ 400, unsteady boundary layer instabilities take place and become progressively dominant.When the driving force is sufficiently large, namely at Ra ≈ 1300 and 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 Eq. ( 28) and Eq. ( 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 Ra = 1255 was used as initial condition for these simulations.As Ra is decreased, the system evolves on two branches (blue lines), both differing from the solution obtained for increasing Ra.
Ra ≈ 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-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 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 Ra, the flow structure (number of rolls) and the transport properties (Sh) differ from those obtained starting from, e.g., a linear temperature distribution or from a solution obtained at lower Ra.As a starting point, Ra = 1255 is used by [53], and Ra is progressively decreased.The system evolves following two distinct branches (Fig. 5a), both differing from the solution obtained for increasing Ra.
Determining the scaling of Sh with Ra in the high-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), [67] and [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 Ra, Sh ∼ Ra.An accurate phenomenological description of the flow and scaling arguments is provided by [14].The linear scaling best fitting the two-dimensional numerical results [54] is and it agrees also with the best known theoretical upper bound, for which Sh ≤ 0.0297Ra [53].The asymptotic scaling proposed by [54] [solid black line in Fig. 5(a)] fits well the numerical results, and it is in agreement with the above mentioned linear predictions: The compensated Sherwood number [Fig.5(b)] 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.5(b)].The best-fitting is in this case provided by [59] 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 Ra → ∞ 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 Ra ≈ 5 × 10 5 , i.e. more than one order of magnitude beyond the threshold identified in two-dimensional flows, and further investigations at Ra ≈ 10 6 are required to confirm this finding.
Resolving the flow equations at the Darcy scale at large 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, [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-wavenumber 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.

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 Sc, defined in (2)], while the density difference produced by the scalar determines the driving force [measured by the Rayleigh number 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 Sh, is discussed here.Fig. 6 Sherwood number (Sh) as a function of Rayleigh-Darcy number (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-Ra scaling [Eq.( 28), grey solid line] are also shown.
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 [31], where the porous medium is modelled as a matrix of aligned squares.They explored different values of porosity (0.36 ≤ ϕ ≤ 0.56) and Schmidt numbers (Sc = 1 and 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 Sh, and the scaling Sh(Ra) appears sublinear when the porosity is increased (ϕ = 0.56).At a low Schmidt numbers (Sc = 1), pore-scale effects on the flow structure, e.g.wavenumber or width of the plumes, are qualitatively similar to those at high Schmidt numbers (Sc = 250).In a complementary study, [70] investigated large Schmidt numbers (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 Ra as where a = 0.011 ± 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 sameRa (achieved with different values of porosity) exhibit very different convective transport properties.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, [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 [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 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 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 Sh.The transition from porous convection to unconfined convection is controlled by two  7(a-c LSC, the velocity in the bulk is non-uniform.Follo from the thermal boundary layers mainly go up an temperature at the cell centre is well mixed with v • T it is confirmed that counter-gradient convectiv the LSC and corner rolls, which is due to the bu between the corner-flow rolls and the LSC (Sugiya 2013).The inclusion of an array of obstacles ha flow structures, as displayed in figure 7(d-i).Con the LSC is suppressed due to the impedance of mixing at the cell centre is less efficient.Thermal p boundary layers can penetrate deep in the bulk, form the regions with strong flows that are formed bet relatively large, the characteristic length scales of than the pore scale, and the flow is chaotic and d as shown in figure 7(d); while when Ra is relativ by large-scale plumes, as shown in figure 7(g).Com is found that the counter-gradient convective heat fl media, which is attributed to the suppression of th dynamics due to the impedance of the obstacle arra 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 Sh with Ra T (not Ra) is reported in Fig. 7(c), where the presence of these two distinct regimes is apparent.For sufficiently large Ra T or high porosity, the classical scaling is recovered (Ra T 1/3 , [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 Sh is proposed: where the Reynolds number Re is computed based on the velocity fluctuations and c = 8 is a fitting parameter.This scaling is proved to be well approximated by Sh ∼ Ra T 0.65 [30]).The transition between these regimes appears clearly in Fig. 7(d), when the compensated Sherwood number (Sh /Ra T −0.3 ) is shown as a function of the boundary layer thickness [δ = H/(2 Sh)] divided by the average pore scale (l s ).The situation is schematically illustrated in Fig. 7(a,b).When the thickness of the thermal boundary layer is comparable to the averaged pore length scale (δ/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].

One-sided convection
The one-sided configuration introduced in Sec.2.3 is representative of natural instances like geological CO 2 sequestration [12] and mixing in groundwater flows [78].In these cases a fluidsaturated porous domain [sketched in Fig. 4(c-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 fingers like structures form [80,81], evolve (i.e., grow and merge) and if the Rayleigh-Darcy number is sufficiently large [Ra > O( 103 )] the system may reach a quasisteady regime.In this phase the dimensionless solute flux 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 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 [82,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 [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 F on the flow parameters during the constant flux regime.

Darcy flows
When the Darcy model is considered [Eq.( 13)-( 15)], the flow is uniquely controlled by the Rayleigh-Darcy number Ra, similarly to the Rayleigh-Bénard case discussed in Sec.3.1.Numerical two-dimensional simulations agree on the value of the flux during the constant flux regime, which was initially determined by [56] to be F = 0.017.(32) This observation has been later confirmed by a number of numerical studies [15-17, 60, 85-87].We refer to [87] for a review of literature scaling laws in the presence of variations to this problem (anisotropy, geochemistry, etc.).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 [15], who performed three-dimensional simulations and estimated the flux to be higher than in the corresponding two-dimensional case (32).These results refer toRa ≤ 9×10 3 , and additional data at larger Rayleigh-Darcy numbers are required to determine at the exact value for 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 rate of scalar dissipation.That period extends, from dimensionless time t ¼ 1 to t ¼ 6, which twice the time that it takes for the fingers to reach m of the domain, indeed highlighting the convecre of the dissolution process.It is interesting that r dissipation rate for the boundary-driven dissoluis approximately twice as large as that for the -fluid model, in analogy with the diffusive flux for nsional diffusion from a one-sided boundary vs two-sided diffusion from an initial sharp uity.mpute the time-averaged mean scalar dissipation uring the time period of constant dissolution flux for the canonical model and t 2 ½2; 8 for gue-fluid model).For simulations with high Ra ), the scalar dissipation rate appears to be indeof Ra (Fig. 2, inset).Given the fluctuations of the lar dissipation rate over time, one cannot reject ypothesis that the dissolution flux is independent yleigh number.Indeed, we fit a power law to hi from the high-resolution simulations as a function is yields a best fit (with 95% confidence bounds) :0120 AE 0:0013ÞRa þ0:031AE0:012 for the canonical nd hi % ð0:0072 AE 0:0012ÞRa À0:017AE0:017 for the -fluid model.These results provide conclusive that the classical Darcy-Boussinesq model of e mixing predicts a regime in which the dissoluand subsequent mixing is constant and, in the high Ra, independent of the Rayleigh number.er, recent experimental studies using analogue ke methanol-ethylene glycol and water (MEG-5] and propylene glycol and water (PG-water) ngering instabilities [23].The scalar dissipay also inform about the time evolution of the th of the flow, which has been shown to exhibit ignature in turbulent flows [24] and heteroges media [25].(4) exposes the fundamental relationship ng rate, dissolution flux, and mean scalar dise, and makes it evident that any Rayleigher-law dependence of the dissolution flux F cted also in the mean scalar dissipation rate hi. to the scalar dissipation rate is particularly aracterize the time evolution in the analogue-In this case, there is no proper dissolution flux a convection-dominated mixing of the two ible fluids.Since all boundaries are no-flux the mean concentration hci is constant, and for the mean square concentration reduces to @ t hc 2 i ¼ À2hi: (5) i provides a fundamentally new way to charmacroscopic evolution of convective mixing, us way to quantify any dependence on Ra. m high-resolution computer simulations of the quations (1), for both the canonical model and e-fluid model.We employ the so-called stream orticity formulation, in which the equation for function and the ADE transport equation are entially at each time step [26].We solve the tion equation using a spectral method [9,27], entration equation using a sixth-order compact ence discretization and a third-order Rungestepping scheme [26].We trigger the densitybility by introducing a small perturbation on ration at the boundary (for the canonical syshorizontal initial interface (for the analogue-), as it is commonly done [9,14,15].ts of a typical simulation are shown in Fig. 1. logy of the convective instability is well known : after an onset period in which a diffusion layer tween the two fluids, the layer develops a oneility in which downward-moving protrusions entially, eventually developing into bloblike thin necks at the roots of the fingers; these interact, merging into each other, and coarsena way that well-developed fingers then attract ed fingers [Fig.1(a)].A snapshot of the simudissipation rate illustrates that the regions ids are actively mixing coincide with the edges ty-driven fingers [Fig.1(b)].This behavior is ualitatively by laboratory experiments with a stem in an Hele-Shaw cell [Fig.1(c)].we plot the time evolution of the mean scalar rate for both the canonical Rayleigh-Be ´nardl and the analogue-fluid model, and for differof Ra.For each case, there is a regime of constant rate of scalar dissipation.That period extends, roughly, from dimensionless time t ¼ 1 to t ¼ 6, which is about twice the time that it takes for the fingers to reach the bottom of the domain, indeed highlighting the convective nature of the dissolution process.It is interesting that the scalar dissipation rate for the boundary-driven dissolution case is approximately twice as large as that for the analogue-fluid model, in analogy with the diffusive flux for one-dimensional diffusion from a one-sided boundary problem vs two-sided diffusion from an initial sharp discontinuity.
We compute the time-averaged mean scalar dissipation rate hi during the time period of constant dissolution flux (t 2 ½1; 6 for the canonical model and t 2 ½2; 8 for the analogue-fluid model).For simulations with high Ra (> 5000), the scalar dissipation rate appears to be independent of Ra (Fig. 2, inset).Given the fluctuations of the mean scalar dissipation rate over time, one cannot reject the null hypothesis that the dissolution flux is independent of the Rayleigh number.Indeed, we fit a power law to hi obtained from the high-resolution simulations as a function of Ra.This yields a best fit (with 95% confidence bounds) hi % ð0:0120 AE 0:0013ÞRa þ0:031AE0:012 for the canonical model, and hi % ð0:0072 AE 0:0012ÞRa À0:017AE0:017 for the analogue-fluid model.These results provide conclusive evidence that the classical Darcy-Boussinesq model of convective mixing predicts a regime in which the dissolution flux and subsequent mixing is constant and, in the range of high Ra, independent of the Rayleigh number.
However, recent experimental studies using analogue fluids, like methanol-ethylene glycol and water (MEGwater) [15] and propylene glycol and water (PG-water) dimension adds significant complexity to the fingering phenomena [15,88], with the flow structure being more complex and dynamical [58].

Pore-scale and Hele-Shaw flows
The determination of F has been carried out beyond the Darcy model via pore-scale simulations and experiments, and via Hele-Shaw setups.
As discussed in Sec.3.1, a classical argument for Darcy convection requires that Sh scales linearly with Ra [67,68].The theoretical interpretation is that in natural convection 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 Ra, i.e., Sh ∼ Ra, it is possible to have an expression for Sh that is independent of H [18,68,90].As a result [see also Eq. ( 24)], the flux F is expected to be independent of Ra, as it emerges from Darcy simulations [e.g., see correlation (32)].Despite this robust theoretical framework, a different scaling for Sh(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 Figs.8(a) and 8(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 ≪ 1), the flow behaves as a laminar Poiseuille flow, i.e., the gapaverage 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 Sh(Ra) scaling was proposed by [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. 8(a)] 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. 8(a)], 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 Sh ∼ 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 Sh to reduce with respect to the classical scaling.An analogue approach was employed by [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. 8(b)].They obtained the scaling Sh ∼ Ra 0.76 and identified the plumes spacing as the key parameter controlling the Sherwood number.Similar results are derived by [91] (Hele-Shaw and beads, Sh ∼ Ra 0.84 ), [92] (Hele-Shaw, Sh ∼ Ra 0.76 ) and [93] (Hele-Shaw, Sh ∼ 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-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 densityconcentration curve, and with a concentrationdependent 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 nonmonotonic density profile is shown in Fig. 8(c).It was found that the dissolution flux is determined by the mean scalar dissipation rate, χ.Mixing in porous media has a universal character, and the non-linear 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) nonmonotonic 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 [90] for non-Boussinesq and compressible flows reveals that the scaling Sh ≈ 181.02 + 0.165Ra 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 non-linear 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 [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 Sh withRa.Later studies [25,96] indicate that within the same value of permeability the scaling Sh ∼ Ra holds.In general, Sh may still be a function of 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 Sh(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 non-linear scaling observed, and will be discussed in detail in the Sec. 5.

Finite-size effects
Domain features like lateral confinement, thickness-induced Hele-Shaw dispersion and porescale dispersion have been identified to play a role on the non-linear scaling of Sh with Ra or the flow structure.The influence of these finite-size effects on convection will be reviewed in this section.

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 flows dynamics.
The flow in a porous Rayleigh-Bénard system at large Ra consists of two distinct regions (see Sec. 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 wavenumber, 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 Ra → ∞ suggests that k ∼ Ra 5/14 , in fair agreement with numerical measurements that give k ∼ Ra 0.4 [54].In three dimensions, theoretical results [99] indicate that k ∼ Ra 1/2 , which is in excellent agreement with numerical measurements of [64] and [58], who obtained k ∼ Ra 0.52 and k ∼ 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 [62] revealed that identifying the wavenumber 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-Ra scaling of the interior inter-plume spacing will require extremely long simulations in very wide computational domains.
In three-dimensions, the effect of the domain confinement has been investigated by [58].They performed numerical simulations at Ra = 10 4 in Rayleigh-Bénard configuration, in domains having variable extension in one of the wall-parallel directions, namely x in Fig. 9(a), 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 A = L/H.Four values of A are considered in Fig. 9, with the domain progressively increasing in size from A = 1/8 to A = 1.The corresponding temperature fields, taken at the centreline (z = 1/2) and close to the bottom wall (z = 0.005), are shown in Figs.9(b)-(i).A strong confinement of the domain presents dramatic effects on the flow structures.For sufficiently large domains, e.g.A = 1, the near-wall cells reported in Fig. 9(b) are randomly oriented and show a wide distribution of sizes.When the domain width is progressively reduced, the cells are strongly constrained [Figs.9(f,h)] and eventually end up in an extremely ordered pattern [Figs.9(d)].The same applies to the flow structures at the centerline that for small domains (A ≤ 1/4) form sheet-like plumes.More quantitative results, estimated by means of the horizontal radial mean wavenumber 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 Ra = 10 4 the flow is independent of the size of the domain for A ≥ 1.
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 [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 threedimensional behaviour.They observed that for sufficiently large Rayleigh-Darcy numbers (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 threedimensional 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 threedimensional 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 co-existence 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. 1(c) relative to geological sequestration of carbon dioxide.Using the model for two-phase gravity currents proposed by [105], [11] analyzed 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 ≥ 140).

Hele-Shaw flows
The working principle of the Hele-Shaw apparatus, briefly introduced in Sec.4.2, is illustrated in Fig. 10(a).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 ≪ 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 ≤ C ≤ C 1 , which produces a maximum density difference ∆ρ within the domain.In this frame, the analogy between Hele-Shaw and Darcy flow has been investigated quantitatively by [97], who observed that a combination of fluid properties (Schmidt number, Sc), cell geometry (anisotropy ratio, ϵ = √ k/H) and flow velocity (U , defined in (8), which depends on 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 valid for ϵ small, Sc ≥ 1 and ϵ 2 Ra ≪ 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 time scale is H/U and the pressure scale is µU H/k.The concentration is made dimensionless as and k is the unit vector with direction opposite to gravity.Eq. ( 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  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 Sc = O( 103 )], 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 ϵ, as predicted by [97], possibly explaining the discrepancy observed between Darcy simulations [21] and Hele-Shaw experiments [20].Finally, the theoretical work proposed by [97] has been recently generalized by [106,107] to the case of more complex systems characterised by the presence of two layers of fluids with nonmonotonic 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 Sh / χ, which is independent of the cell geometry (ϵ) and directly proportional to Ra.Using this theoretical framework, they suggest that a possible reason for the sublinear scaling observed by [20] is the flow regime (Hele-Shaw regime) in which the experiments are performed.

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 Sec.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 Sec. 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 [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 work, most of 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 [109], magnetic resonance imaging (MRI) of threedimensional convective flows in opaque media were presented, and plumes at low Rayleigh-Darcy numbers (< 20π 2 ) were visualized.Also X-ray computed tomography (CT) imaging scan is now frequently used to study mixing of miscible fluids.[110,111] provided correlations for Sherwood as a function of Péclet and Rayleigh-Darcy number, and observed a sublinear scaling for Sh with Ra, with exponent 0.40 and 0.93 respectively.The same methodology was employed by [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 Sh ∼ Ra holds, in contrast with previous studies.This discrepancy may be due to the relatively short range and small values of 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, [114] achieved larger Rayleigh-Darcy numbers (≤ 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(Ra) holds.
These works agree upon the fact that dispersion is crucial in determining the Sh(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, [87,117] observed that the flow structures are influenced by the strength of dispersion and the dissolution rate 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, [26] proposed that, in addition to the Rayleigh-Darcy number, a flow with dispersion is controlled by a dispersive Rayleigh-Darcy number with D 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 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 L /D T (see Sec. 2.2.1) is also important to determine the flow character.As a result, the parameter space for convective porous media flows with dispersion is controlled by at least three parameters: Ra, Ra d and r.In order to quantify the relative importance of molecular diffusion to transverse dispersion, one can introduce the parameter [84,118] that can be used to rewrite the dispersion tensor (16) in dimensionless form as: This expression suggests that the case of pure diffusion is recovered when D T ≪ D, corresponding to ∆ ≫ 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 α L = D L /U ≈ d, we can rewrite Eq. (35) as In bead packs, the permeability can be inferred from the Kozeny-Carman correlation [45,119], i.e.
where k C = 5 is the Carman constant for monodispersed spheres randomly packed [120].As a result, we can provide an expression for Ra d and Ra that is an explicit function of the domain (g, H), fluid (D, µ, ∆ρ) and medium (ϕ, d, r) properties.This information is particularly important when we characterise the flow in the three-dimensional parameters space (Ra d ,Ra, r), which we will do in the following.First, we reduce the parameters space to (Ra,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 Sh occurs for r > 10, provided that Ra and Ra d are sufficiently large (namely, 10 4 and 10 3 , respectively [84]).For r ≤ 1, the  flow is qualitatively similar or that observed in the absence of dispersion [54].With respect to the remaining parameters, Ra and Ra d , we can rewrite both as a function of the beads diameter and find thatRa d ∼ 1/d andRa ∼ 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 (Ra,Ra d ) shown in Fig. 11, corresponding to Ra d ∼ Ra −1/2 .Using realistic laboratory properties, we obtain that a possible range for the experimental parameters (Ra d ,Ra) at variable d consists of the circles in Fig. 11 (each series of circles corresponds to one value of density difference, ∆ρ).Alternatively, we consider the case in which the medium is fixed (d constant) and the fluid density contrast varies.Since Ra d is independent of any fluid property, a variation of ∆ρ 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 ∆.It follows that this is achieved when (∆ρ) −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 ∆ρ and d.
With the aid of numerical simulations the problem of decoupling two of the governing flow parameters (Ra,Ra d ) can be solved, and their relative effect on the Sh or F can be investigated.[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 ∆, which determines the mechanism controlling convection.If ∆ > 10 5 , the flow is ruled by molecular diffusion, plumes grow symmetrically [Fig.12(a)] and the structure is analogue to the symmetric flow observed in the absence of dispersion [see Fig. 4(b-i)].When ∆ < 1, mechanical dispersion dominates over convection, and its inherent anisotropy (r ≫ 1) sets the non-symmetric flow structure (fan flow) shown in Fig. 12(c), in which plumes widen as they move away from the wall.A similar behaviour is observed in the corresponding one-sided cases [Fig.12(b,d)] by [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 non-zero velocity parallel to the wall is admitted (free-slip), which produces solute spread due to dispersion.Alternatively, 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 (Ra,Ra d ) space, where regions controlled by different mechanisms are separated by dashed lines.[26] found that in Rayleigh-Bénard configuration, when ∆ > O(1) molecular diffusion dominates over mechanical dispersion, although a small contribution of mechanical dispersion may increase Sh.When 0.02 < ∆ < O(1), both mechanical dispersion and molecular diffusion determine the value of Sh.A linear scaling Sh ∼ Ra holds when ∆ < 0.02, but Ra d is also important since it determines the prefactor of the scaling law, as it has been later observed by [121].
These results for dispersion-dominated flows (∆ < 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 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, ∆ sets the flow behaviour: diffusion dominated [∆ < O( 1 believe that such one-to-one comparison between results of dispersion for Hele-Shaw and bead packs flows may not be appropriate. Recent works investigated the role of mechanical dispersion with the aid of simulations.A possible complementary approach with respect to the formulation used by [26,84] is proposed by [123], where the dispersive Rayleigh number is replaced by a parameter quantifying the strength of longitudinal dispersion compared to molecular diffusion.More recently, [118] performed simulations in one-sided configuration.They modelled fluids with constant viscosity and linear densityconcentration profiles, and derived a linear correlation between Sh and Ra, where the prefactor is a function of molecular to dispersive Rayleigh-Darcy numbers,Ra /Ra d .This correlation fits well their results, but it does not fully capture the trend predicted by [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 (Ra,Ra d ) space are required.
Finally, an novel approach consists of including the effects of dispersion also in the momentum equation.[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 √ 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 porescale effects in advection-dominated systems in the absence of convection [125].

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 non-linear phenomenon taking place in multiphase and multiscale systems, eventually located thousands of meters 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 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 Sh ∼ Ra.While in the two-dimensional case [54,62] this scaling sets in at Ra ≤ 10 4 , in three-dimensional flows [58,59] the ultimate state is expected to take place at Ra ≥ 5 × 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 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 F is quasi steady.While Darcy simulations report a constant Ra-independent value F [16,17,56,60], experiments in bead packs [19] and Hele-Shaw cells [20] revealed that F is a function of 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 (Ra), (ii) convection and dispersion (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 Sec.6.1 we will briefly review recent advancements in experimental techniques, and in Sec.6.2 we will also discuss the importance of additional effects not considered in previous sections of this paper.

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 [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 nonintrusive 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 µm and in time of 6 s.Recent experiments [131,132] have shown that the resolution can be further lowered down to 2.3 µ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 [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 lab models is essential to capture the role of the medium heterogeneities on the solute mixing.The technique proposed by [139] addresses this issue, and it consists of a cell made of 3Dprinted elementary blocks designed to be easily re-arranged 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 easy 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.

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 Ra ρr /∆ρ ≫ 1, being ∆ρ the maximum density difference and ρr the reference fluid density.Interestingly, this condition is independent of ∆ρ and it is widely fulfilled for geothermal processes, when Ra ≈ 10 1 − 10 3 and ρ r /∆ρ ≈ 10 2 −10 3 [147].Numerical simulations of the fullycompressible 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 highand 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 nontrivial manners.[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.[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].
The liquid phase in the secondary mushy layer is further enriched in component C and depleted in component B.Even for the relatively simple phase diagram shown in figure 4, one can also envision adjusting the initial liquid compo-

Fig. 1
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).Saltrich (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.

Fig. 2
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 (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.

Fig. 4
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 (u • n = 0), and side boundaries may be also considered as periodic.The reference frame (x, z) and gravity (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 twodimensional simulations at 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 χ, F and dt⟨C 2 ⟩/2 is reported (the operator dt stands for the time-derivative).Quantities are computed as in Eq. (25) and made dimensionless with respect to the length-scale L = ℓ.The time-averaged value of F is also shown (dashed lines) in panels (b-ii) and (c-ii).

Fig. 5
Fig. 5 Darcy simulations.(panel a) Sherwood number (Sh) as a function of Rayleigh-Darcy number (Ra), and (panel b) in compensated form (Sh /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 Eq. (28) and Eq.(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 Ra = 1255 was used as initial condition for these simulations.As Ra is decreased, the system evolves on two branches (blue lines), both differing from the solution obtained for increasing Ra.

FIGURE 7 .Fig. 7
FIGURE 7. Typical snaps velocity magnitude |u| and direction at (a-c) (Ra, ) (10 7 , 0.82).Circles in (d temperature equation (2.2), Snapshot of the concentration c at dimensionless 1 from a simulation of the analogue-fluid system with 00 and constant viscosity (R ¼ 0).A computational 2 Â 1536 cells was used, and only a small window of ation domain is shown.(b) Corresponding snapshot of dissipation rate , for the same simulation as that in , dark color corresponds to high values of , and regions of active mixing.(c) Snapshot of a surrogate alar dissipation rate ¼ rc Á D m rc (obtained from sity) from a laboratory experiment with a PG-water a Hele-Shaw cell, illustrating that mixing is primarily to narrow layers along the edges of the horizontal and the density-driven fingers.

FIG. 1 .Fig. 8
FIG. 1.(a) Snapshot of the concentration c at dimensionless time t ¼ 1 from a simulation of the analogue-fluid system with Ra ¼ 10000 and constant viscosity (R ¼ 0).A computational grid of 512 Â 1536 cells was used, and only a small window of the simulation domain is shown.(b) Corresponding snapshot of the scalar dissipation rate , for the same simulation as that in (a).Here, dark color corresponds to high values of , and indicates regions of active mixing.(c) Snapshot of a surrogate of the scalar dissipation rate ¼ rc Á D m rc (obtained from light intensity) from a laboratory experiment with a PG-water system in a Hele-Shaw cell, illustrating that mixing is primarily confined to narrow layers along the edges of the horizontal interface and the density-driven fingers.

Fig. 9
Fig. 9 Influence of lateral confinement (domain width) on the development of the flow structures.Three-dimensional Rayleigh-Bénard simulations performed at Ra = 10 4 are shown [adapted with permission from 58, 59].(a) Dimensionless temperature distribution (ϑ) in a cubic domain, being 0 and 1 the values at top and bottom boundaries, respectively, with gravity 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 A = L/H increases from A = 1/8 (panels b,c) to A = 1 (panels h,i).For each value of A, 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 colorbars apply to centreline and near-wall panels.
[97]: (i) Darcy regime [Fig.10(b)] when ϵ → 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.10(c)] when ϵ ≪ 1, ϵ 2 Ra ≪ 1 and Sc ≥ 1, characterised a gradient of concentration across the cell gap, but with one single finger; and (iii) three-dimensional regime [Fig.10(d)], 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,[97] provided an evidence for the reduction of the scaling exponent (discussed in Sec.4.2) for convective flows in the Hele-Shaw regime.These finding were later confirmed by the laboratory experiments of[25], where the flux has

H
Fig.10(a) Front view of 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 (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[97] are shown.

Fig. 11
Fig. 11 Parameters space (Ra,Ra d ) with indication of iso-∆ lines (blue lines).With respect to experiments in bead packs, if only d varies, the parameters (Ra,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 (∆ρ).Conversely, if only ∆ρ is varied in the experiments, (Ra,Ra d ) are locked to horizontal lines (diamonds).Triangles indicate the slope of the green and blue lines.

FIG. 1 .
FIG. 1. Time-averaged horizontal-mean concentration profile C and snapshots of the concentration field C from DNS at Ram = 20 000 and r = 10 for different Rad .The domain aspect ratio is L = 5.In panel (a), only half of C is shown due to its antisymmetry about the midplane, and the z values on the horizontal axis are nonuniformly spaced to clearly show the structure near the wall.Increasing mechanical disperison (decreasing Rad ) thickens the diffusive boundary layer, coarsens the flow pattern, and stabilizes the flow.Moreover, the convection transitions to a fan-flow structure at Rad < 5000.so that Rad becomes the only parameter controlling the dynamics of the system for fixed r.Thus, at large Ram the concentration field C and the buoyancy velocity u are determined solely by the dispersive Rayleigh number Rad , as confirmed by our DNS data.Once C and u become invariant in the limit of Ram → ∞, Fm ∼ c1, and Fd ∼ c2Ram with the constants c1 and c2 determined by Rad , as shown in Fig. 4(a).123801-6
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, ∆ sets the flow behaviour: diffusion dominated [∆ < O(1)], dispersion dominated [∆ > 0.02] or influenced by both diffusion and dispersion [0.02 < ∆ < O(1)].