Hostname: page-component-cd9895bd7-7cvxr Total loading time: 0 Render date: 2024-12-23T01:13:51.256Z Has data issue: false hasContentIssue false

Towards the understanding of convective dissolution in confined porous media: thin bead pack experiments, two-dimensional direct numerical simulations and physical models

Published online by Cambridge University Press:  16 May 2024

Marco De Paoli*
Affiliation:
Physics of Fluids Group and Max Planck Center for Complex Fluid Dynamics and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands Institute of Fluid Mechanics and Heat Transfer, TU Wien, 1060 Vienna, Austria
Christopher J. Howland*
Affiliation:
Physics of Fluids Group and Max Planck Center for Complex Fluid Dynamics and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands
Roberto Verzicco*
Affiliation:
Physics of Fluids Group and Max Planck Center for Complex Fluid Dynamics and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands Dipartimento di Ingegneria Industriale, University of Rome ‘Tor Vergata’, 00133 Roma, Italy Gran Sasso Science Institute, 67100 L'Aquila, Italy
Detlef Lohse*
Affiliation:
Physics of Fluids Group and Max Planck Center for Complex Fluid Dynamics and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands Max Planck Institute for Dynamics and Self-Organization, 37077 Gottingen, Germany

Abstract

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

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2024. Published by Cambridge University Press.

1. Introduction

Natural convective flows take place when density differences within a fluid domain drive the fluid motion. Natural convection is the dominant heat and mass transport mechanism for many flows of practical interest, both in nature (Cardin & Olson Reference Cardin and Olson1994; Marshall & Schott Reference Marshall and Schott1999; Hartmann, Moy & Fu Reference Hartmann, Moy and Fu2001) and in industrial applications (Krepper, Hicken & Jaegers Reference Krepper, Hicken and Jaegers2002; Bejan Reference Bejan2013). One of the particularly relevant ones in geophysics is the case of convection in porous media (De Paoli Reference De Paoli2023), in which the fluid occupies the interstitial space of a solid porous matrix. This configuration is typical of subsurface transport phenomena, and is relevant for petroleum migration (Simmons, Fenstemaker & Sharp Reference Simmons, Fenstemaker and Sharp2001), water contamination (LeBlanc Reference LeBlanc1984; Van Der Molen & Van Ommen Reference Van Der Molen and Van Ommen1988), underground hydrogen storage (Zivar, Kumar & Foroozesh Reference Zivar, Kumar and Foroozesh2021; Krevor et al. Reference Krevor, de Coninck, Gasda, Ghaleigh, de Gooyert, Hajibeygi, Juanes, Neufeld, Roberts and Swennenhuis2023), superficial formations in dry salty lakes (Lasser, Ernst & Goehring Reference Lasser, Ernst and Goehring2021; Lasser et al. Reference Lasser, Nield, Ernst, Karius, Wiggs, Threadgold, Beaume and Goehring2023) and sea ice growth (Wettlaufer, Worster & Huppert Reference Wettlaufer, Worster and Huppert1997; Feltham et al. Reference Feltham, Untersteiner, Wettlaufer and Worster2006; Middleton et al. Reference Middleton, Thomas, de Wit and Tison2016), to name a few.

Recently, this problem has been the subject of extended investigations due to the implications it can bear for geological carbon sequestration (Huppert & Neufeld Reference Huppert and Neufeld2014; Emami-Meybodi et al. Reference Emami-Meybodi, Hassanzadeh, Green and Ennis-King2015). In this process, carbon dioxide (CO$_{2}$) is injected into underground formations located 1–3 km beneath the Earth's surface, with the aim of permanent storage. After injection, carbon dioxide remains in a supercritical state due to the high pressure at the injections depths, but its density is lower compared with that of the resident fluid (brine), and the injected volume of carbon dioxide migrates on top of the brine layer (MacMinn et al. Reference MacMinn, Neufeld, Hesse and Huppert2012). This situation is undesired, because it may cause a further migration of carbon dioxide to the upper layers, eventually leading to a return of CO$_{2}$ into the atmosphere. However, carbon dioxide is partially soluble in brine, and when these fluids mix a heavier solution (CO$_{2}$ + brine) forms. When a sufficiently thick layer of this denser mixture builds up, it may become unstable and form convective instabilities (De Paoli Reference De Paoli2021). These flow structures have a favourable effect on the efficiency of the CO$_{2}$ dissolution mechanism, since they contribute to fluid mixing (i.e. CO$_{2}$ permanently stored in brine) in a much more effective manner compared with pure diffusion (Ennis-King, Preston & Paterson Reference Ennis-King, Preston and Paterson2005; Xu, Chen & Zhang Reference Xu, Chen and Zhang2006). On the other hand, the process is made complex by the presence of these nonlinear instabilities, and making reliable predictions of the dissolution dynamics is a non-trivial task.

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

The mixing process of CO$_{2}$ and brine is idealised considering an homogeneous porous slab with a constant concentration difference applied at the horizontal boundaries (Rayleigh–Bénard), or as a semi-infinite (one-sided) domain where the concentration is fixed at top and no flux is applied at bottom boundary (Hewitt Reference Hewitt2020). In the Rayleigh–Bénard–Darcy case, the system attains a statistically steady state, which has been predicted theoretically (Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948) and recently accurately described numerically in two (Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2012; De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2016; Wen et al. Reference Wen, Akhbari, Zhang and Hesse2018; Ulloa & Letelier Reference Ulloa and Letelier2022) and three dimensions (Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2014; De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022a; Dhar et al. Reference Dhar, Meunier, Nadal and Méheust2022), also at high Rayleigh–Darcy numbers, $\operatorname {\mathit {Ra}}^*$ (the Rayleigh–Darcy number indicates the relative importance of convective and diffusive contributions in a flow through a porous medium). The one-sided-Darcy configuration is characterised by a transient behaviour (Riaz et al. Reference Riaz, Hesse, Tchelepi and Orr2006; Fu, Cueto-Felgueroso & Juanes Reference Fu, Cueto-Felgueroso and Juanes2013) consisting of three main phases: (i) an initial diffusive regime in which the mixing layer grows and becomes eventually unstable (De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2017), (ii) a convective phase characterised by a constant dissolution rate (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012) and (iii) a finite-size regime when the strength of convection reduces gradually (Slim et al. Reference Slim, Bandi, Miller and Mahadevan2013).

The Rayleigh–Bénard–Darcy and the one-sided-Darcy systems have been well characterised numerically and theoretically in case of thermal convection at the Darcy scale, i.e. when the flow structures are large compared with the scale of the pores (Hewitt Reference Hewitt2020). Numerical results suggest that a linear relationship exists between the dimensionless mass transfer coefficient (Sherwood number, $\operatorname {\mathit {Sh}}$) and the Rayleigh–Darcy number ($\operatorname {\mathit {Ra}}^*$), in both the Rayleigh–Bénard–Darcy case (Hewitt et al. Reference Hewitt, Neufeld and Lister2012; Pirozzoli et al. Reference Pirozzoli, De Paoli, Zonta and Soldati2021) and during the convective regime of the one-sided-Darcy system (Slim Reference Slim2014; De Paoli et al. Reference De Paoli, Zonta and Soldati2017). Experimental measurements in these configurations (Neufeld et al. Reference Neufeld, Hesse, Riaz, Hallworth, Tchelepi and Huppert2010; Backhaus, Turitsyn & Ecke Reference Backhaus, Turitsyn and Ecke2011; De Paoli, Alipour & Soldati Reference De Paoli, Alipour and Soldati2020; Brouzet, Méheust & Meunier Reference Brouzet, Méheust and Meunier2022), however, suggest a different qualitative and also quantitative behaviour compared with the corresponding Darcy simulations, with a nonlinear scaling of $\operatorname {\mathit {Sh}}$ with $\operatorname {\mathit {Ra}}^*$. This discrepancy is likely due to non-Darcy effects (Liang et al. Reference Liang, Wen, Hesse and DiCarlo2018), i.e. to the pore-induced flow dynamics not captured by Darcy simulations. Therefore, resolving the flow and the solute transport at the pore scale is crucial to make reliable models to incorporate in large-scale simulations and to predict the underground fluids migration and mixing, and it represents the main motivation for this work.

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

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

A porous Rayleigh–Taylor system consists of two miscible fluids of different density initially arranged in an unstable configuration and immersed in a porous matrix. After an initial diffusive phase (Wang et al. Reference Wang, Nakanishi, Hyodo and Suekane2016, Reference Wang, Nakanishi, Teston and Suekane2018), the flow is driven by convection, and due to the transient nature of this system, the mixing rate of the two fluids varies in time. In geophysical applications, it is crucial to determine the mixing rate of the involved fluids, to be able to make reliable predictions on the evolution of the volume of subsurface flow. The porous Rayleigh–Taylor system has been studied at the Darcy scale with the aid of Hele-Shaw cells and Darcy simulations (De Wit Reference De Wit2020). The flow behaviour is quantified by the mixing length, i.e. the vertical extension of the mixing region. The growth rate of the mixing length is controlled by the combined action of buoyancy (density difference) and drag (viscous dissipation through the medium) (Boffetta & Musacchio Reference Boffetta and Musacchio2022; De Paoli et al. Reference De Paoli, Perissutti, Marchioli and Soldati2022b). In the case of miscible fluids, the role of diffusion across the fluid–fluid interface is also crucial (Gopalakrishnan et al. Reference Gopalakrishnan, Carballido-Landeira, De Wit and Knaepen2017), as it weakens convection by reducing local density gradients.

As a result of this time-dependent finger growth process, the mixing dynamics is also transient, and it may be quantified via the mean scalar dissipation rate (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012). The mean scalar dissipation is particularly convenient because it can be exactly related to other flow quantities, e.g. the averaged concentration fluctuations, and it has been already described numerically at the Darcy scale (De Paoli et al. Reference De Paoli, Giurgiu, Zonta and Soldati2019a; De Paoli Reference De Paoli2023). However, in the Rayleigh–Taylor case, the behaviour of dispersion cannot be captured by Darcy simulations despite its crucial role, as it can influence dramatically the onset of the gravitational instabilities (Menand & Woods Reference Menand and Woods2005) and the mixing dynamics. In this flow configuration, the full flow dynamics has been studied with the aid of three-dimensional simulations by Sardina et al. (Reference Sardina, Brandt, Boffetta and Mazzino2018), where the authors consider a thermal convection at relatively large values of porosity ($0.6\unicode{x2013}1$). They proposed a model to incorporate the effect of the medium within a friction coefficient to be included in the Navier–Stokes equations. The case of solutal convection in geological formations may be markedly different, since the porosity is low (0.2–0.4) and the Schmidt numbers about two orders of magnitude larger than in the thermal case, and we aim precisely at this gap.

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

Our findings are relevant to the convection of miscible fluids in porous media. However, we note that some differences occur when CO$_2$ sequestration is considered, including the non-monotonic relationship between the fluid density and the solute concentration (Hidalgo et al. Reference Hidalgo, Dentz, Cabeza and Carrera2015), and the partial miscibility of the phases involved (Huppert & Neufeld Reference Huppert and Neufeld2014). These effects, as well as additional limitations (dimensionality of the systems and idealised structure of the media) later discussed in detail, prevent a direct application of the present findings to CO$_2$ dissolution in geological formations. Nonetheless, these results represent a fundamental ingredient required to build a modelling framework for large-scale simulations, as they are obtained in a well-defined and controlled physical system.

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

2. Problem formulation

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

(2.1a,b)\begin{equation} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}=0,\quad \frac{\partial C}{\partial n}=0 ,\end{equation}

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

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

In the present flow configuration, the dimensionless parameters controlling the system can be grouped into three main categories: medium parameters (Darcy number, porosity), fluid parameters (Schmidt number) and flow parameters (Rayleigh, Rayleigh–Darcy, Péclet and Reynolds numbers).

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

(2.2)\begin{equation} \operatorname{\mathit{Da}}=\frac{k}{H^2},\end{equation}

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

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

(2.3)\begin{equation} \operatorname{\mathit{Sc}}=\frac{\mu}{\rho_0 D}.\end{equation}

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

(2.4)\begin{equation} U=\frac{g\Delta\rho k}{\mu},\end{equation}

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

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

(2.5)\begin{equation} \ell=\frac{\phi D}{U}.\end{equation}

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

(2.6a,b)\begin{equation} \operatorname{\mathit{Ra}}=\frac{g\Delta\rho H^3}{\mu D},\quad \operatorname{\mathit{Ra}}_d=\frac{g\Delta\rho d^3}{\mu D}=\operatorname{\mathit{Ra}}\left(\frac{d}{H}\right)^3, \end{equation}

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

(2.7)\begin{equation} \operatorname{\mathit{Ra}}^*=\frac{g\Delta\rho k H}{\mu\phi D}=\frac{UH}{\phi D}=\frac{H}{\ell}=\frac{\operatorname{\mathit{Ra}}\operatorname{\mathit{Da}}}{\phi}.\end{equation}

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

(2.8a,b)\begin{equation} Re=\frac{\operatorname{\mathit{Ra}}^*\operatorname{\mathit{Da}}^{1/2}}{Sc}\ll 1 ,\quad \operatorname{\mathit{Pe}}=\operatorname{\mathit{Ra}}^*\operatorname{\mathit{Da}}^{1/2}\ll 1,\end{equation}

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

The flow scales of the experiments are listed in table 1, whereas the dimensionless parameters corresponding to present experiments and simulations are reported in table 2.

Table 1. Dimensional parameters of the performed experiments. The domain size is constant ($L=H=200$ mm, $B=28.5$ mm). Fluid temperature ($\vartheta$), solute mass fraction ($\omega _M$) and density contrast ($\Delta \rho$) are reported. Fluid viscosity and diffusivity are $\mu =9.2\times 10^{-4}\ {\rm Pa}\ {\rm s}^{-1}$ and $D=1.65\times 10^{-9}\ {\rm m}^2\ {\rm s}^{-1}$, respectively. Water density is assumed to be $\rho _0=10^3\ {\rm kg}\ {\rm m}^{-3}$. Bead diameter ($d$), medium porosity ($\phi$) and permeability ($k$) are indicated. Flow length scale $\ell$ and velocity scale $U$, defined as in (2.5) and (2.4), respectively, are also reported.

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

3. Methodology

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

3.1. Laboratory experiments

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

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

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

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

To allow here a reliable comparison of experimental results against the direct numerical simulations, the fluids employed in the experiments have been specifically selected because of the linear dependency of the density of the solution with respect to the solute concentration (Slim et al. Reference Slim, Bandi, Miller and Mahadevan2013; De Paoli et al. Reference De Paoli, Perissutti, Marchioli and Soldati2022b). This condition is indeed met not only in the simulations presented here, but also in most of the numerical works available in the literature.

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

(3.1)\begin{equation} \omega(C) = \frac{C}{\rho(C)}.\end{equation}

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

(3.2)\begin{equation} \rho=\rho_0\left[1+\frac{\Delta\rho }{\rho_0C_0}(C-C_0)\right].\end{equation}

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

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

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

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

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

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

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

Figure 3. Processing of images. (a) Normalised intensity profile $I_1$, defined in (3.4), shown over the entire measurement region for the experiment E12. The initial position of the interface ($z/H=0$, red solid line) is reported. The instantaneous position of the interface, defined as an isocontour of $I_{1}$, is also shown (blue line). (b) Corrected intensity profile $I_2$, defined as in (3.5). The corrected intensity $I_{2}$ is lower than 1 within the mixing region. (c) Horizontally averaged intensity profiles. The lower edge of the mixing region, determined by the position where $\bar {I}_{2}$ is lower than a given threshold value, is reported (dashed red line). This is related to the mixing length $h$, discussed in detail in § 5. See also supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.328.

The preprocessed fields are analysed to produce intensity profiles and infer the time-dependent evolution of the interface. The preprocessing consists of several steps, illustrated in figure 3. The initial light intensity distribution $I(x,z,t=0)$, obtained from the raw images, is used to compute and store the initial mean light intensity values of the high ($I_H$, $z/H>0$) and low ($I_L$, $z/H<0$) fluid density layers. The normalised light intensity field

(3.4)\begin{equation} I_1(x,z,t) = \frac{I(x,z,t)-I_H}{I_L-I_H},\end{equation}

is suitable to visualise the instantaneous flow configuration. It is used, for instance, to identify the instantaneous position of the interface of the mixing region (figure 3a). However, we observe in figure 3(c) that the horizontally averaged intensity profile $\bar {I}_1$ varies smoothly within the mixing region, and therefore it is not a good indicator to determine the edge of the interface. We introduce the corrected intensity (figure 3b), $I_2$, defined as

(3.5)\begin{equation} I_2(x,z,t) = 2[\tfrac{1}{2} + I_1(x,z,t) - \hat{I}_1(x,z,0) ],\end{equation}

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

3.2. Numerical simulations

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

(3.6)$$\begin{gather} \partial_t \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u} ={-}\rho_0^{{-}1}\boldsymbol{\nabla} p + \nu \nabla^2 \boldsymbol{u} - g\beta C \hat{\boldsymbol{z}}, \end{gather}$$
(3.7)$$\begin{gather}\partial_t C + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})C = D \nabla^2 C, \end{gather}$$

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

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

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

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

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

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

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

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

(3.8a,b)\begin{equation} C(z)|_{t=0} = C_0[\mathcal{H}(z) + A(z) W] , \quad A(z) = 0.1 \mathrm{sech}^2(100 z/H), \end{equation}

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

4. Flow dynamics

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

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

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

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

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

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

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

Figure 7. Snapshots of the concentration field from simulations S10, S12 and S6. In the first two simulations, $H/d=70$ is fixed, and only the Rayleigh number $Ra_d$ is varied from $10^5$ to $10^6$ (left to centre). In the third simulation (S6, right column), $Ra_d=10^6$ matches the centre column but $H/d$ is decreased to 35, so the obstacles are larger compared with the domain height. Each snapshot is presented at a given fraction of the time scale $H/U$, with time increasing top to bottom (as specified in the left column). Videos showing the evolution of the concentration field in selected simulations (S4, S6 and S12) are available as supplementary movies (movie 2, movie 3 and movie 4, respectively).

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

5. Mixing length and concentration profiles

The mixing length is defined as the vertical extent of the region where solute concentration is non-uniform. Multiple ways to quantify the mixing length have been proposed (Cook, Cabot & Miller Reference Cook, Cabot and Miller2004), depending on whether a bulk-focused measure is taken or it tracks the spread of the fastest growing fingers (e.g. $0.01\le \bar {C}/C_0 \le 0.99$; see Gopalakrishnan et al. Reference Gopalakrishnan, Carballido-Landeira, De Wit and Knaepen2017). The latter method is used here to determine the mixing length obtained from the experiments, where a threshold value corresponding to $0.96$ is applied to the horizontally averaged corrected light intensity as defined in (3.5), $\bar {I}_2$. The mixing length is then computed assuming that the flow is symmetric with respect to the domain centreline, and a time correction is also applied (see § A.3).

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

(5.1)\begin{equation} \frac{\bar{C}(z,t)}{C_0} = \begin{cases} 0 & z\le-h/2 ,\\ 1/2 + z/h & |z| < h/2, \\ 1 & z\ge h/2. \end{cases} \end{equation}

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

(5.2)\begin{equation} h = \frac{6}{C_0^2}\int_{-\infty}^{\infty} \bar{C} (1 - \bar{C}) \,\mathrm{d}z ,\end{equation}

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

The dimensionless mixing length scaled with respect to $\ell$ is shown for experiments (colourbar) and simulations (legend) in figure 8. Simulations are presented for both the regular bead pattern (solid lines) and the perturbed bead pattern (dotted lines) as shown in figure 4. Asymptotically, both experiments and simulations follow the scaling $h=2Ut$ (dashed line). The buoyancy velocity $U$ represents the terminal velocity of a rising (falling) parcel of light (heavy) fluid surrounded by heavy (light) fluid, and it is achieved at the equilibrium between the driving force, represented by buoyancy, and the drag due to viscous forces within the medium. This model is a simplified representation of the flow, and any diffusion effect or interaction with the flow structures is neglected. However, it provides a good first-order estimate for the growth of the mixing region. The observed behaviour gives a precise indication of the flow dynamics: after a starting phase in which the flow is influenced by the initial condition, the late-stage dynamics is controlled by convection, with the mixing layer growing at the characteristic buoyancy velocity $U$.

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

In three-dimensional porous systems, Sardina et al. (Reference Sardina, Brandt, Boffetta and Mazzino2018) showed that the mixing layer grows linearly for low values of porosity ($\phi =0.6$), whereas it follows the classical turbulent quadratic scaling when the porosity approaches 1. This suggests that a linear growth of the mixing region is expected also for the values of porosity considered in this study. However, the additional degree of freedom provided by the third spatial dimension may have an effect on the velocity at which the plumes spread: compared with the two-dimensional case, more pathways will be available to the individual fingers, which can spread more horizontally. As a result, we expect that during the convective regime the mixing length will still grow at a rate that is constant in time and lower compared with the two-dimensional case. A uniform growth of the mixing layer is also observed in two- and three-dimensional Darcy flow at large Rayleigh–Darcy numbers (Boffetta, Borgnino & Musacchio Reference Boffetta, Borgnino and Musacchio2020; Borgnino, Boffetta & Musacchio Reference Borgnino, Boffetta and Musacchio2021). Also in this case, when comparing the evolution of two- and three-dimensional flows some differences emerge. The growth rate of the mixing length is larger in two dimensions than in three dimensions (Borgnino et al. Reference Borgnino, Boffetta and Musacchio2021), and the transition between these regime occurs sharply (Boffetta & Musacchio Reference Boffetta and Musacchio2022), as soon as the domain thickness is larger than the wavelength of the most unstable mode. The nature of this transition in pore-resolved flows has not been explored yet. At the Darcy scale and at moderate Rayleigh–Darcy numbers, a superlinear scaling has been observed (De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2019b), which has been interpreted has a finite size effect (Boffetta & Musacchio Reference Boffetta and Musacchio2022; De Paoli et al. Reference De Paoli, Perissutti, Marchioli and Soldati2022b), i.e. the time/space available for the fingers before touching the horizontal boundaries is insufficient to reach this asymptotic regime.

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

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

6. Flow structure and wavenumber

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

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

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

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

(6.1)\begin{equation} \lambda(t) = \frac{\displaystyle \int_0^\infty k^{{-}1} |\hat{C}(k,t)|^2 \,\mathrm{d}k}{\displaystyle \int_0^\infty |\hat{C}(k,t)|^2 \,\mathrm{d}k}.\end{equation}

As shown in figure 11, the finger structures at the centre of the domain exhibit a coarsening behaviour, with the mean wavelength increasing over time. This coarsening is consistent across all the simulations at varying $Ra_d$ and $H/d$. After initial transient phase related to the onset of the Rayleigh–Taylor instability, the wavelength exhibits a power-law scaling of $t^{1/2}$. Such a scaling is often associated with a diffusive process, although in figure 11 the collapse is observed in terms of the advective scale $U$ rather than the diffusivity $D$. This result of a scaling close to $t^{1/2}$ is in fair agreement with the Darcy simulations of De Paoli et al. (Reference De Paoli, Giurgiu, Zonta and Soldati2019a), where the mean wavenumber $\kappa \sim 1/\lambda$ exhibits a $t^{-0.6}$ scaling during the growth of the mixing layer. Boffetta et al. (Reference Boffetta, Borgnino and Musacchio2020) also observe $t^{1/2}$ scaling for the horizontal wavenumber in three-dimensional Darcy simulations of Rayleigh–Taylor instability, finding a collapse with $\lambda \sim \sqrt {Dt}$. This diffusive behaviour is attributed to nonlinear processes, including plume merging and coarsening, and it is not a simple direct consequence of molecular diffusion. This contrasts to our pore-scale simulations, where the rate of coarsening is independent of the molecular diffusivity. As the mixing layer grows in our simulations, and the fingers coarsen to scales larger than a few multiples of the pore scale, the Darcy assumption seems to become more relevant, as would be expected.

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

7. Mean scalar dissipation

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

(7.1)\begin{equation} \chi = D \langle |\boldsymbol{\nabla} C|^2 \rangle, \end{equation}

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

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

The mean scalar dissipation can be related to other global quantities of the flow (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; De Paoli Reference De Paoli2023), such as the dissolution rate (in case of systems permeable to solute at the boundaries) or the volume-averaged squared concentration, $\langle C^2\rangle$. Given our no flux boundary conditions for concentration, the relation

(7.2)\begin{equation} \partial_t \langle C^2\rangle ={-}2 \chi, \end{equation}

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

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

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

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

(7.3)\begin{equation} C = \frac{C_0}{2}\left[1 +\mathrm{erf} \left( \frac{z}{2 \sqrt{D t}}\right)\right] . \end{equation}

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

(7.4)\begin{equation} \chi = D \langle |\boldsymbol{\nabla} C|^2 \rangle = \frac{ D}{H} \int_{{-}H/2}^{H/2} |\partial_{z} C|^2\, \mathrm{d}z \approx \frac{ D}{H} \int_{-\infty}^{\infty} |\partial_{z} C|^2\, \mathrm{d}z = \sqrt{\frac{ D}{8 {\rm \pi}t}} \frac{C_0^{2}}{H}.\end{equation}

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

(7.5)\begin{equation} \frac{\phi\chi}{UC_0^{2}/H} =\frac{1}{\sqrt{8{\rm \pi}}} \left(\frac{t}{\phi\ell/U}\right)^{{-}1/2}.\end{equation}

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

(7.6)\begin{equation} \chi = \sqrt{\frac{ D}{8 {\rm \pi}t}} \frac{C_0^{2}}{H}g(\xi), \quad\text{with } g(\xi) =\text{erf}(\xi)-\text{erf} \left(\frac{\xi}{\sqrt{2}}\right)\sqrt{2\, {\rm e}^{-\xi^2}},\end{equation}

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

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

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

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

(7.7)\begin{equation} \chi = D \langle |\boldsymbol{\nabla} C |^2 \rangle \approx D \frac{h}{H} \langle |\boldsymbol{\nabla} C|^2 \rangle_{ML} , \end{equation}

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

(7.8)\begin{equation} |\boldsymbol{\nabla} C| \approx \frac{C_0}{2\sqrt{{\rm \pi} D t}} . \end{equation}

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

(7.9)\begin{equation} \chi \approx D \frac{2 U t}{H} \frac{C_0^{2}}{4{\rm \pi} D t} = \frac{1}{2{\rm \pi}} \frac{UC_0^{2}}{H} . \end{equation}

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

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

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

(7.10)\begin{equation} \frac{\bar{C}(z,t)}{C_0} = \begin{cases} 0 & \text{for } -H/2\le z\le -h/2,\\ 1/2+z/(2Ut) & \text{for } |z|< h/2,\\ 1 & \text{for } h/2\le z\le H/2, \end{cases}\end{equation}

which is nothing but (5.1) with $h=2Ut$. Using (7.2) with approximation (7.10), we estimate the mean scalar dissipation as

(7.11)\begin{equation} \chi \approx{-}\frac{1}{2} \frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{{-}H/2}^{{-}h(t)/2}\bar{C}^2 \,\mathrm{d}z+\int_{{-}h(t)/2}^{h(t)/2}\bar{C}^2 \,\mathrm{d}z+\int_{h(t)/2}^{H/2}\bar{C}^2 \,\mathrm{d}z\right) = \frac{1}{6} \frac{UC_0^{2}}{H} . \end{equation}

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

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

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

(7.12)\begin{equation} L_{i}(t) =L+2\frac{L}{\lambda(t)}h(t).\end{equation}

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

(7.13)\begin{equation} \chi = D \langle |\boldsymbol{\nabla} C |^2 \rangle = \frac{DL_{i}}{HL}\int_{-\delta/2}^{+\delta/2}|\partial_{n}C|^{2}\mathop{}\!\mathrm{d} n \approx \frac{D L_i}{HL} \frac{C_0^2 \delta}{4{\rm \pi} Dt}. \end{equation}

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

(7.14)\begin{equation} \delta = \beta \sqrt{dU t},\end{equation}

with $\beta$ being a fitting parameter. Note that the effective diffusivity $dU$ determined here coincides with the effective longitudinal dispersion usually considered for bead packs (Liang et al. Reference Liang, Wen, Hesse and DiCarlo2018; Tsinober, Shavit & Rosenzweig Reference Tsinober, Shavit and Rosenzweig2023). We now substitute the expressions obtained in (7.12) and (7.14) into the definition (7.13), together with the evolution of the mixing length $h=2Ut$ and the wavelength measured from the simulations $\lambda /d=\alpha \sqrt {t/(d/U)}$, as in figure 11. Finally, we obtain an expression for $\chi (t)$ in the convective regime:

(7.15)\begin{equation} \chi=\frac{C_0^{2}U}{H}\frac{\beta}{\alpha{\rm \pi}} \left[1+\frac{\alpha}{4}\left(\frac{t}{d/U}\right)^{{-}1/2}\right],\end{equation}

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

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

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

We observe that the behaviour of $\chi$ shown in figure 14 is remarkably dissimilar compared with the corresponding Darcy simulations (De Paoli et al. Reference De Paoli, Giurgiu, Zonta and Soldati2019a), where it is observed to increase with time up to impingement of the fingers on the horizontal boundaries. Hidalgo et al. (Reference Hidalgo, Dentz, Cabeza and Carrera2015) observed that in presence of fluids characterised by non-monotonic density–concentration profiles, which leads to a considerably different flow topology, the scalar dissipation remains constant during the convective regime. We speculate that the pore-induced dispersion effects (Dentz, Icardi & Hidalgo Reference Dentz, Icardi and Hidalgo2018) are likely responsible for these differences.

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

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

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

8. Summary, conclusions and outlook

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

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

In this study we have focused on modelling the mixing dynamics and the flow structure of porous media convection. We used simulations performed in two-dimensional arrays of circular objects and experiments consisting of thin three-dimensional packs of spherical beads (Hele-Shaw type geometry). In a porous Rayleigh–Taylor system, at the Darcy scale, the mixing region is observed to grow faster in two-dimensional domains than in three-dimensional ones (Borgnino et al. Reference Borgnino, Boffetta and Musacchio2021). Unlike in the turbulent case, the transition between these regimes occurs sharply when the thickness of the domain exceeds the wavelength of the most unstable mode (Boffetta & Musacchio Reference Boffetta and Musacchio2022). The nature of this transition in pore-resolved flows has not been explored yet, and in the future it will be of interest to extend the present study to simulations in three dimensions, as well as two-dimensional experiments (De et al. Reference De, Singh, Fulcrand, Méheust, Meunier and Nadal2022), to allow for a one-to-one comparison of these findings and to investigate the effect of the dimensionality of the flow on the evolution of a pore-resolved system. At the pore scale, three-dimensionality provides more pathways for fingers to percolate through, and the interfacial area of three-dimensional fingers will be greater than for two-dimensional fingers. Such effects may potentially allow for greater dispersion, and quantifying the associated impact of the solute on the larger-scale spread will be a key focus of future simulations.

The present findings are relative to domains having a dimension up to few hundred pores. In contrast, geophysical applications involve formations that may be orders of magnitude larger (Huppert & Neufeld Reference Huppert and Neufeld2014), and modelling of the entire dissolution process may be required (Szulczewski, Hesse & Juanes Reference Szulczewski, Hesse and Juanes2013; Wang, Vuik & Hajibeygi Reference Wang, Vuik and Hajibeygi2022). For instance, in the case of carbon sequestration it is desired to determine the time required to dissolved a prescribed fraction of the injected fluid volume, or to estimate the horizontal spread of the current of CO$_2$ (MacMinn et al. Reference MacMinn, Neufeld, Hesse and Huppert2012; De Paoli Reference De Paoli2021). Performing pore-scale studies of such flows is beyond the present capabilities and Darcy simulations including the effect of dispersion are required (Dentz, Hidalgo & Lester Reference Dentz, Hidalgo and Lester2023). In this context, one remarkable finding of the present study is the behaviour of the mean scalar dissipation during the convection-dominated regime, which is observed to decrease in time in pore-resolved flows (see figure 14) whereas it grows in Darcy flows (De Paoli et al. Reference De Paoli, Giurgiu, Zonta and Soldati2019a). A different behaviour is observed by Hidalgo et al. (Reference Hidalgo, Dentz, Cabeza and Carrera2015) in the presence of fluids characterised by non-monotonic density–concentration profiles. In their Darcy simulations, this property of the fluids leads to a remarkably different flow topology with a nearly flat interface between the two fluid layers, and the scalar dissipation remains constant during the convective regime. We speculate that the pore-induced dispersion effects (Dentz et al. Reference Dentz, Icardi and Hidalgo2018), not captured by Darcy simulations, are likely responsible for these differences and should be accounted by suitable dispersion models developed from pore-resolved data.

The results presented are relevant to homogeneous and isotropic porous media at the porosity $\phi =0.37$. Geological formations are generally characterised by a more complex structure, typically heterogeneous at multiple scales (Hidalgo, Neuweiler & Dentz Reference Hidalgo, Neuweiler and Dentz2022), anisotropic (Ennis-King et al. Reference Ennis-King, Preston and Paterson2005) and with porosity values that can be as low as 0.05 (Bickle et al. Reference Bickle, Kampman, Chapman, Ballentine, Dubacq, Galy, Sirikitputtisak, Warr, Wigley and Zhou2017). Exploring the effect of such irregularities of the medium on the flow structure would represent an important extension of this study to be considered for future works.

Supplementary movies

Supplementary movies are available at https://doi.org/10.1017/jfm.2024.328.

Funding

This project has received funding from the European Union's Horizon Europe research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 101062123. The work of CJH was funded by the Max Planck Center for Complex Fluid Dynamics. We acknowledge PRACE for awarding us access to MareNostrum4 at the Barcelona Supercomputing Center (BSC), Spain and IRENE at Très Grand Centre de Calcul (TGCC) du CEA, France (project 2021250115). This work was also carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. This research was funded in part by the Austrian Science Fund (FWF) [grant J-4612].

Declaration of interests

The authors report no conflict of interest.

Data availability statement

The data presented in this work are openly available at Howland & De Paoli (Reference Howland and De Paoli2024).

Appendix A. Additional experimental details

A.1. Fluid density

Following the empirical correlations of Novotný & Söhnel (Reference Novotný and Söhnel1988), the density of an aqueous solution of KMnO$_4$ can be determined as a function of the solute concentration $C$ (expressed in mol dm$^{-3}$) and the temperature of the solution $\vartheta$ (expressed in $^{\circ }\text {C}$) as

(A1)\begin{equation} \rho(C,\vartheta)=\rho_w(\vartheta)+ A_1(\vartheta)C+A_2(\vartheta)C^{3/2}, \end{equation}

where the water density $\rho _w(\vartheta )$ and the temperature-dependent coefficients $A_1(\vartheta )$ and $A_2(\vartheta )$, are given by

(A2)\begin{align} \rho_w(\vartheta) &={+}9.997\times 10^2\frac{\text{kg}}{\text{m}^3}+2.044\times 10^{{-}1}\frac{\text{kg}}{^\circ\text{C} \ \text{m}^3 }\cdot \vartheta\nonumber\\ &\quad -6.174\times 10^{{-}2}\frac{\text{kg}}{(^\circ\text{C})^{3/2} \ \text{m}^3 }\cdot \vartheta^{3/2}, \end{align}
(A3)\begin{align}A_1(\vartheta) &={+}1.223\times 10^{{-}1}\frac{\text{kg}}{\text{mol}}-1.029\times 10^{{-}4}\frac{\text{kg}}{^\circ\text{C}\text{ mol}}\cdot \vartheta\nonumber\\ &\quad +8.093\times 10^{{-}6}\frac{\text{kg}}{(^\circ\text{C})^2\text{ mol}}\cdot \vartheta^2, \end{align}
(A4)\begin{align} A_2(\vartheta) &={-}1.485\times 10^{1}\frac{\text{kg}}{\text{m}^3}\left(\frac{\text{dm}^3}{\text{mol}}\right)^{3/2}+9.079\times 10^{{-}1}\frac{\text{kg}}{^\circ\text{C }\text{m}^3}\left(\frac{\text{dm}^3}{\text{mol}}\right)^{3/2}\cdot \vartheta\nonumber\\&\quad -7.566\times 10^{{-}3}\frac{\text{kg}}{(^\circ\text{C }^2)\text{m}^3}\left(\frac{\text{dm}^3}{\text{mol}}\right)^{3/2}\cdot \vartheta^2. \end{align}

The relative dependence of density of the solution, $\rho$, mass fraction, $\omega$, and solute concentration, $C$, is shown in figure 17 (blue solid lines). We report here $C(\omega )$ (figure 17a), $\rho (\omega )-\rho _w$ (figure 17b) and $\rho (C)-\rho _w$ (figure 17c), being $\rho _w$ the water density defined as in (A2). The density of the mixture (figure 17c), $\rho (C)$, is well approximated by a linear function of the solute concentration (3.2) (dashed grey line).

Figure 17. Given the definition (3.1) and the empirical correlations (A1)–(A4), the relative dependence of density of the solution, $\rho$, mass fraction, $\omega$, and solute concentration, $C$, is obtained (blue solid lines). We report $C(\omega )$ (a), $\rho (\omega )-\rho _w$ (b) and $\rho (C)-\rho _w$ (c), where $\rho _w$ is the water density defined as in (A2). The profiles shown here correspond to $\vartheta =25\,^{\circ }\text {C}$. The fluid properties are well approximated by linear functions (grey dashed lines).

A.2. Characterisation of the porous medium

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

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

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

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

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

References

Ataei-Dadavi, I., Rounaghi, N., Chakkingal, M., Kenjeres, S., Kleijn, C.R. & Tummers, M.J. 2019 An experimental study of flow and heat transfer in a differentially side heated cavity filled with coarse porous media. Intl J. Heat Mass Transfer 143, 118591.CrossRefGoogle Scholar
Backhaus, S., Turitsyn, K. & Ecke, R.E. 2011 Convective instability and mass transport of diffusion layers in a Hele-Shaw geometry. Phys. Rev. Lett. 106 (10), 104501.CrossRefGoogle Scholar
Bejan, A. 2013 Convection Heat Transfer. John Wiley & Sons.CrossRefGoogle Scholar
Bickle, M., Kampman, N., Chapman, H., Ballentine, C., Dubacq, B., Galy, A., Sirikitputtisak, T., Warr, O., Wigley, M. & Zhou, Z. 2017 Rapid reactions between CO$_2$, brine and silicate minerals during geological carbon storage: modelling based on a field CO$_2$ injection experiment. Chem. Geol. 468, 1731.CrossRefGoogle Scholar
Boffetta, G., Borgnino, M. & Musacchio, S. 2020 Scaling of Rayleigh–Taylor mixing in porous media. Phys. Rev. Fluids 5 (6), 62501.CrossRefGoogle Scholar
Boffetta, G. & Mazzino, A. 2017 Incompressible Rayleigh–Taylor turbulence. Annu. Rev. Fluid Mech. 49, 119143.CrossRefGoogle Scholar
Boffetta, G. & Musacchio, S. 2022 Dimensional effects in Rayleigh–Taylor mixing. Phil. Trans. R. Soc. 380 (2219), 20210084.CrossRefGoogle ScholarPubMed
Borgnino, M., Boffetta, G. & Musacchio, S. 2021 Dimensional transition in Darcy–Rayleigh–Taylor mixing. Phys. Rev. Fluids 6 (7), 074501.CrossRefGoogle Scholar
Brouzet, C., Méheust, Y. & Meunier, P. 2022 Co 2 convective dissolution in a three-dimensional granular porous medium: an experimental study. Phys. Rev. Fluids 7 (3), 033802.CrossRefGoogle Scholar
Cardin, P. & Olson, P. 1994 Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core. Phys. Earth Planet. Inter. 82 (3-4), 235259.CrossRefGoogle Scholar
Chakkingal, M., Kenjereš, S., Ataei-Dadavi, I., Tummers, M.J. & Kleijn, C.R. 2019 Numerical analysis of natural convection with conjugate heat transfer in coarse-grained porous media. Intl J. Heat Fluid Flow 77, 4860.CrossRefGoogle Scholar
Cook, A.W., Cabot, W. & Miller, P.L. 2004 The mixing transition in Rayleigh–Taylor instability. J. Fluid Mech. 511, 333362.CrossRefGoogle Scholar
De, N., Singh, N., Fulcrand, R., Méheust, Y., Meunier, P. & Nadal, F. 2022 Two-dimensional micromodels for studying the convective dissolution of carbon dioxide in 2D water-saturated porous media. Lab on a Chip 22 (23), 46454655.CrossRefGoogle ScholarPubMed
De Paoli, M. 2021 Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide. Phys. Fluids 33 (1), 016602.CrossRefGoogle Scholar
De Paoli, M. 2023 Convective mixing in porous media: a review of Darcy, pore-scale and Hele-Shaw studies. Eur. Phys. J. E 46 (12), 129.CrossRefGoogle ScholarPubMed
De Paoli, M., Alipour, M. & Soldati, A. 2020 How non-Darcy effects influence scaling laws in Hele-Shaw convection experiments. J. Fluid Mech. 892, A41.CrossRefGoogle Scholar
De Paoli, M., Giurgiu, V., Zonta, F. & Soldati, A. 2019 a Universal behavior of scalar dissipation rate in confined porous media. Phys. Rev. Fluids 4 (10), 101501.CrossRefGoogle Scholar
De Paoli, M., Perissutti, D., Marchioli, C. & Soldati, A. 2022 b Experimental assessment of mixing layer scaling laws in Rayleigh–Taylor instability. Phys. Rev. Fluids 7 (9), 093503.CrossRefGoogle Scholar
De Paoli, M., Pirozzoli, S., Zonta, F. & Soldati, A. 2022 a Strong Rayleigh–Darcy convection regime in three-dimensional porous media. J. Fluid Mech. 943, A51.CrossRefGoogle Scholar
De Paoli, M., Zonta, F. & Soldati, A. 2016 Influence of anisotropic permeability on convection in porous media: implications for geological CO$_2$ sequestration. Phys. Fluids 28 (5), 056601.CrossRefGoogle Scholar
De Paoli, M., Zonta, F. & Soldati, A. 2017 Dissolution in anisotropic porous media: modelling convection regimes from onset to shutdown. Phys. Fluids 29 (2), 026601.CrossRefGoogle Scholar
De Paoli, M., Zonta, F. & Soldati, A. 2019 b Rayleigh–Taylor convective dissolution in confined porous media. Phys. Rev. Fluids 4, 023502.CrossRefGoogle Scholar
De Wit, A. 2020 Chemo-hydrodynamic patterns and instabilities. Annu. Rev. Fluid Mech. 52, 531555.CrossRefGoogle Scholar
Delgado, J.M.P.Q. 2007 Longitudinal and transverse dispersion in porous media. Chem. Engng Res. Des. 85 (9), 12451252.CrossRefGoogle Scholar
Dentz, M., Hidalgo, J.J. & Lester, D. 2023 Mixing in porous media: concepts and approaches across scales. Transp. Porous Med. 146 (1–2), 553.CrossRefGoogle Scholar
Dentz, M., Icardi, M. & Hidalgo, J.J. 2018 Mechanisms of dispersion in a porous medium. J. Fluid Mech. 841, 851882.CrossRefGoogle Scholar
Dhar, J., Meunier, P., Nadal, F. & Méheust, Y. 2022 Convective dissolution of carbon dioxide in two-and three-dimensional porous media: the impact of hydrodynamic dispersion. Phys. Fluids 34 (6), 064114.CrossRefGoogle Scholar
Dullien, F.A.L. 1991 Porous Media: Fluid Transport and Pore Structure. Academic.Google Scholar
Eckel, A.M., Liyanage, R., Kurotori, T. & Pini, R. 2022 Spatial moment analysis of convective mixing in three-dimensional porous media using X-ray CT images. Ind. Engng Chem. Res. 62 (1), 762–774.Google Scholar
Emami-Meybodi, H., Hassanzadeh, H., Green, C.P. & Ennis-King, J. 2015 Convective dissolution of CO$_2$ in saline aquifers: progress in modeling and experiments. Intl J. Greenh. Gas Control 40, 238266.CrossRefGoogle Scholar
Ennis-King, J., Preston, I. & Paterson, L. 2005 Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys. Fluids 17 (8), 084107.CrossRefGoogle Scholar
Fadlun, E.A., Verzicco, R., Orlandi, P. & Mohd-Yusof, J. 2000 Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. J. Comput. Phys. 161 (1), 3560.CrossRefGoogle Scholar
Feltham, D.L., Untersteiner, N., Wettlaufer, J.S. & Worster, M.G. 2006 Sea ice is a mushy layer. Geophys. Res. Lett. 33 (14).CrossRefGoogle Scholar
Fu, X., Cueto-Felgueroso, L. & Juanes, R. 2013 Pattern formation and coarsening dynamics in three-dimensional convective mixing in porous media. Phil. Trans. R. Soc. A 371 (2004), 20120355.CrossRefGoogle ScholarPubMed
Gasow, S., Kuznetsov, A.V., Avila, M. & Jin, Y. 2021 A macroscopic two-length-scale model for natural convection in porous media driven by a species-concentration gradient. J. Fluid Mech. 926, A8.CrossRefGoogle Scholar
Gasow, S., Kuznetsov, A.V. & Jin, Y. 2022 Prediction of pore-scale-property dependent natural convection in porous media at high Rayleigh numbers. Intl J. Therm. Sci. 179, 107635.CrossRefGoogle Scholar
Gasow, S., Lin, Z., Zhang, H.C., Kuznetsov, A.V., Avila, M. & Jin, Y. 2020 Effects of pore scale on the macroscopic properties of natural convection in porous media. J. Fluid Mech. 891, A25.CrossRefGoogle Scholar
Gopalakrishnan, S.S., Carballido-Landeira, J., De Wit, A. & Knaepen, B. 2017 Relative role of convective and diffusive mixing in the miscible Rayleigh–Taylor instability in porous media. Phys. Rev. Fluids 2, 012501.CrossRefGoogle Scholar
Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: a unifying theory. J. Fluid Mech. 407, 2756.CrossRefGoogle Scholar
Grossmann, S. & Lohse, D. 2001 Thermal convection for large Prandtl numbers. Phys. Rev. Lett. 86 (15), 3316.CrossRefGoogle ScholarPubMed
Happel, J. & Brenner, H. 2012 Low Reynolds Number Hydrodynamics: With Special Applications to Particulate Media, vol. 1. Springer Science & Business Media.Google Scholar
Hartmann, D.L., Moy, L.A. & Fu, Q. 2001 Tropical convection and the energy balance at the top of the atmosphere. J. Clim. 14 (24), 44954511.2.0.CO;2>CrossRefGoogle Scholar
Haughey, D.P. & Beveridge, G.S.G. 1969 Structural properties of packed beds - a review. Can. J. Chem. Engng 47 (2), 130140.CrossRefGoogle Scholar
Hewitt, D.R. 2020 Vigorous convection in porous media. Proc. Math. Phys. Engng Sci. 476 (2239), 20200111.Google ScholarPubMed
Hewitt, D.R., Neufeld, J.A. & Lister, J.R. 2012 Ultimate regime of high Rayleigh number convection in a porous medium. Phys. Rev. Lett. 108 (22), 224503.CrossRefGoogle Scholar
Hewitt, D.R., Neufeld, J.A. & Lister, J.R. 2014 High Rayleigh number convection in a three-dimensional porous medium. J. Fluid Mech. 748, 879895.CrossRefGoogle Scholar
Hidalgo, J.J., Dentz, M., Cabeza, Y. & Carrera, J. 2015 Dissolution patterns and mixing dynamics in unstable reactive flow. Geophys. Res. Lett. 42 (15), 63576364.CrossRefGoogle Scholar
Hidalgo, J.J., Fe, J., Cueto-Felgueroso, L. & Juanes, R. 2012 Scaling of convective mixing in porous media. Phys. Rev. Lett. 109 (26), 264503.CrossRefGoogle ScholarPubMed
Hidalgo, J.J., Neuweiler, I. & Dentz, M. 2022 Advective trapping in the flow through composite heterogeneous porous media. Transp. Porous Med. 143 (3), 599618.CrossRefGoogle Scholar
Horton, C.W. & Rogers, F.T. Jr. 1945 Convection currents in a porous medium. J. Appl. Phys. 16 (6), 367370.CrossRefGoogle Scholar
Howland, C.J. & De Paoli, M. 2024 Data supporting ‘Towards the understanding of convective dissolution in confined porous media: thin bead pack experiments, two-dimensional direct numerical simulations and physical models’. Available at: https://doi.org/10.4121/897ba0bb-c3e5-4e31-9e6a-31f6a19f2e6c.CrossRefGoogle Scholar
Huppert, H.E. & Neufeld, J.A. 2014 The fluid mechanics of carbon dioxide sequestration. Annu. Rev. Fluid Mech. 46, 255272.CrossRefGoogle Scholar
Korba, D. & Li, L. 2022 Effects of pore scale and conjugate heat transfer on thermal convection in porous media. J. Fluid Mech. 944, A28.CrossRefGoogle Scholar
Krepper, E., Hicken, E.-F. & Jaegers, H. 2002 Investigation of natural convection in large pools. Intl J. Heat Fluid Flow 23 (3), 359365.CrossRefGoogle Scholar
Krevor, S., de Coninck, H., Gasda, S.E., Ghaleigh, N.S., de Gooyert, V., Hajibeygi, H., Juanes, R., Neufeld, J., Roberts, J.J. & Swennenhuis, F. 2023 Subsurface carbon dioxide and hydrogen storage for a sustainable energy future. Nat. Rev. Earth Environ. 4 (2), 102118.CrossRefGoogle Scholar
Lapwood, E.R. 1948 Convection of a fluid in a porous medium. Math. Proc. Camb. Phil. Soc. 44, 508521.CrossRefGoogle Scholar
Lasser, J., Ernst, M. & Goehring, L. 2021 Stability and dynamics of convection in dry salt lakes. J. Fluid Mech. 917, A14.CrossRefGoogle Scholar
Lasser, J., Nield, J.M., Ernst, M., Karius, V., Wiggs, G.F.S., Threadgold, M.R., Beaume, C. & Goehring, L. 2023 Salt polygons and porous media convection. Phys. Rev. X 13 (1), 011025.Google Scholar
LeBlanc, D.R. 1984 Sewage Plume in a Sand and Gravel Aquifer, Cape Cod, Massachusetts, vol. 2218. US Geological Survey.Google Scholar
Liang, Y., Wen, B., Hesse, M.A. & DiCarlo, D. 2018 Effect of dispersion on solutal convection in porous media. Geophys. Res. Lett. 45 (18), 96909698.CrossRefGoogle Scholar
Liu, S., Jiang, L., Chong, K.L., Zhu, X., Wan, Z.H., Verzicco, R., Stevens, R.J.A.M., Lohse, D. & Sun, C. 2020 From Rayleigh–Bénard convection to porous-media convection: how porosity affects heat transfer and flow structure. J. Fluid Mech. 895, A18.CrossRefGoogle Scholar
MacMinn, C.W., Neufeld, J.A., Hesse, M.A. & Huppert, H.E. 2012 Spreading and convective dissolution of carbon dioxide in vertically confined, horizontal aquifers. Water Resour. Res. 48 (11).CrossRefGoogle Scholar
Marshall, J. & Schott, F. 1999 Open-ocean convection: observations, theory, and models. Rev. Geophys. 37 (1), 164.CrossRefGoogle Scholar
Menand, T. & Woods, A.W. 2005 Dispersion, scale, and time dependence of mixing zones under gravitationally stable and unstable displacements in porous media. Water Resourc. Res. 41 (5).CrossRefGoogle Scholar
Middleton, C.A., Thomas, C., de Wit, A. & Tison, J.-L. 2016 Visualizing brine channel development and convective processes during artificial sea-ice growth using Schlieren optical methods. J. Glaciol. 62 (231), 117.CrossRefGoogle Scholar
Neufeld, J.A., Hesse, M.A., Riaz, A., Hallworth, M.A., Tchelepi, H.A. & Huppert, H.E. 2010 Convective dissolution of carbon dioxide in saline aquifers. Geophys. Res. Lett. 37 (22).CrossRefGoogle Scholar
Nield, D. & Bejan, A. 2017 Convection in Porous Media, 5th edn. Springer.CrossRefGoogle Scholar
Novotný, P. & Söhnel, O. 1988 Densities of binary aqueous solutions of 306 inorganic substances. J. Chem. Engng Data 33 (1), 4955.CrossRefGoogle Scholar
Ostilla-Monico, R., Yang, Y., van der Poel, E.P., Lohse, D. & Verzicco, R. 2015 A multiple-resolution strategy for direct numerical simulation of scalar turbulence. J. Comput. Phys. 301, 308321.CrossRefGoogle Scholar
Pirozzoli, S., De Paoli, M., Zonta, F. & Soldati, A. 2021 Towards the ultimate regime in Rayleigh–Darcy convection. J. Fluid Mech. 911, R4.CrossRefGoogle Scholar
van der Poel, E.P., Ostilla-Mónico, R., Donners, J. & Verzicco, R. 2015 A pencil distributed finite difference code for strongly turbulent wall-bounded flows. Comput. Fluids 116, 1016.CrossRefGoogle Scholar
Riaz, A., Hesse, M., Tchelepi, H.A. & Orr, F.M. 2006 Onset of convection in a gravitationally unstable diffusive boundary layer in porous media. J. Fluid Mech. 548, 87111.CrossRefGoogle Scholar
Sardina, G., Brandt, L., Boffetta, G. & Mazzino, A. 2018 Buoyancy-driven flow through a bed of solid particles produces a new form of Rayleigh–Taylor turbulence. Phys. Rev. Lett. 121 (22), 224501.CrossRefGoogle Scholar
Simmons, C.T., Fenstemaker, T.R. & Sharp, J.M. Jr. 2001 Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges. J. Contam. Hydrol. 52 (1–4), 245275.CrossRefGoogle ScholarPubMed
Slim, A.C. 2014 Solutal-convection regimes in a two-dimensional porous medium. J. Fluid Mech. 741, 461491.CrossRefGoogle Scholar
Slim, A.C., Bandi, M.M., Miller, J.C. & Mahadevan, L. 2013 Dissolution-driven convection in a Hele–Shaw cell. Phys. Fluids 25 (2), 024101.CrossRefGoogle Scholar
Szulczewski, M.L., Hesse, M.A. & Juanes, R. 2013 Carbon dioxide dissolution in structural and stratigraphic traps. J. Fluid Mech. 736, 287315.CrossRefGoogle Scholar
Tsinober, A., Shavit, U. & Rosenzweig, R. 2023 Numerical investigation of the influence of hydrodynamic dispersion on solutal natural convection. Water Resourc. Res. 59 (5), e2023WR034475.CrossRefGoogle Scholar
Ulloa, H.N. & Letelier, J.A. 2022 Energetics and mixing of thermally driven flows in Hele-Shaw cells. J. Fluid Mech. 930, A16.CrossRefGoogle Scholar
Van Der Molen, W.H. & Van Ommen, H.C. 1988 Transport of solutes in soils and aquifers. J. Hydrol. 100 (1), 433451.CrossRefGoogle Scholar
Verzicco, R. 2023 Immersed boundary methods: historical perspective and future outlook. Annu. Rev. Fluid Mech. 55 (1), 129155.CrossRefGoogle Scholar
Verzicco, R. & Orlandi, P. 1996 A finite-difference scheme for three-dimensional incompressible flows in cylindrical coordinates. J. Comput. Phys. 123 (2), 402414.CrossRefGoogle Scholar
Wang, L., Nakanishi, Y., Hyodo, A. & Suekane, T. 2016 Three-dimensional structure of natural convection in a porous medium: effect of dispersion on finger structure. Intl J. Greenh. Gas Control 53, 274283.CrossRefGoogle Scholar
Wang, L., Nakanishi, Y., Teston, A.D. & Suekane, T. 2018 Effect of diffusing layer thickness on the density-driven natural convection of miscible fluids in porous media: modeling of mass transport. J. Fluid Sci. Technol. 13 (1), JFST0002.CrossRefGoogle Scholar
Wang, Y., Vuik, C. & Hajibeygi, H. 2022 Analysis of hydrodynamic trapping interactions during full-cycle injection and migration of CO$_2$ in deep saline aquifers. Adv. Water Resour. 159, 104073.CrossRefGoogle Scholar
Wen, B., Akhbari, D., Zhang, L. & Hesse, M.A. 2018 Convective carbon dioxide dissolution in a closed porous medium at low pressure. J. Fluid Mech. 854, 5687.CrossRefGoogle Scholar
Wettlaufer, J.S., Worster, M.G. & Huppert, H.E. 1997 The phase evolution of Young Sea Ice. Geophys. Res. Lett. 24 (10), 12511254.CrossRefGoogle Scholar
Woods, A.W. 2015 Flow in Porous Rocks. Cambridge University Press.Google Scholar
Xu, P. & Yu, B. 2008 Developing a new form of permeability and Kozeny–Carman constant for homogeneous porous media by means of fractal geometry. Adv. Water Resour. 31 (1), 7481.CrossRefGoogle Scholar
Xu, X., Chen, S. & Zhang, D. 2006 Convective stability analysis of the long-term storage of carbon dioxide in deep saline aquifers. Adv. Water Resour. 29 (3), 397407.CrossRefGoogle Scholar
Zaman, E. & Jalali, P. 2010 On hydraulic permeability of random packs of monodisperse spheres: direct flow simulations versus correlations. Phys. A: Stat. Mech. Appl. 389 (2), 205214.CrossRefGoogle Scholar
Zhong, J., Liu, S. & Sun, C. 2023 On the thermal effect of porous material in porous media Rayleigh–Bénard convection. Flow 3, E13.CrossRefGoogle Scholar
Zivar, D., Kumar, S. & Foroozesh, J. 2021 Underground hydrogen storage: a comprehensive review. Intl J. Hydrog. Energy 46 (45), 2343623462.CrossRefGoogle Scholar
Figure 0

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

Figure 1

Table 1. Dimensional parameters of the performed experiments. The domain size is constant ($L=H=200$ mm, $B=28.5$ mm). Fluid temperature ($\vartheta$), solute mass fraction ($\omega _M$) and density contrast ($\Delta \rho$) are reported. Fluid viscosity and diffusivity are $\mu =9.2\times 10^{-4}\ {\rm Pa}\ {\rm s}^{-1}$ and $D=1.65\times 10^{-9}\ {\rm m}^2\ {\rm s}^{-1}$, respectively. Water density is assumed to be $\rho _0=10^3\ {\rm kg}\ {\rm m}^{-3}$. Bead diameter ($d$), medium porosity ($\phi$) and permeability ($k$) are indicated. Flow length scale $\ell$ and velocity scale $U$, defined as in (2.5) and (2.4), respectively, are also reported.

Figure 2

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

Figure 3

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

Figure 4

Figure 3. Processing of images. (a) Normalised intensity profile $I_1$, defined in (3.4), shown over the entire measurement region for the experiment E12. The initial position of the interface ($z/H=0$, red solid line) is reported. The instantaneous position of the interface, defined as an isocontour of $I_{1}$, is also shown (blue line). (b) Corrected intensity profile $I_2$, defined as in (3.5). The corrected intensity $I_{2}$ is lower than 1 within the mixing region. (c) Horizontally averaged intensity profiles. The lower edge of the mixing region, determined by the position where $\bar {I}_{2}$ is lower than a given threshold value, is reported (dashed red line). This is related to the mixing length $h$, discussed in detail in § 5. See also supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.328.

Figure 5

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

Figure 6

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

Figure 7

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

Figure 8

Figure 7. Snapshots of the concentration field from simulations S10, S12 and S6. In the first two simulations, $H/d=70$ is fixed, and only the Rayleigh number $Ra_d$ is varied from $10^5$ to $10^6$ (left to centre). In the third simulation (S6, right column), $Ra_d=10^6$ matches the centre column but $H/d$ is decreased to 35, so the obstacles are larger compared with the domain height. Each snapshot is presented at a given fraction of the time scale $H/U$, with time increasing top to bottom (as specified in the left column). Videos showing the evolution of the concentration field in selected simulations (S4, S6 and S12) are available as supplementary movies (movie 2, movie 3 and movie 4, respectively).

Figure 9

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

Figure 10

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

Figure 11

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

Figure 12

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

Figure 13

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

Figure 14

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

Figure 15

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

Figure 16

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

Figure 17

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

Figure 18

Figure 17. Given the definition (3.1) and the empirical correlations (A1)–(A4), the relative dependence of density of the solution, $\rho$, mass fraction, $\omega$, and solute concentration, $C$, is obtained (blue solid lines). We report $C(\omega )$ (a), $\rho (\omega )-\rho _w$ (b) and $\rho (C)-\rho _w$ (c), where $\rho _w$ is the water density defined as in (A2). The profiles shown here correspond to $\vartheta =25\,^{\circ }\text {C}$. The fluid properties are well approximated by linear functions (grey dashed lines).

Figure 19

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

Figure 20

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

Supplementary material: File

De Paoli et al. supplementary movie 1

(a) Normalised intensity profile, I1, shown over the entire measurement region for the experiment E12. The initial position of the interface (z/H=0, red solid line) is reported. The instantaneous position of the interface, defined as an iso-contour of I1, is also shown. The color of the iso-contour varies with time. In addition to the current iso-contour, few instantaneous iso-contours at different time instants are reported. (b) Corrected intensity profile I2, which is lower than 1 within the mixing region. (c) Horizontally-averaged intensity profiles. The lower edge of the mixing region, setting the extension of the mixing length h and determined by the position where Ī_2 is lower than a given threshold value, is reported (dashed red line).
Download De Paoli et al. supplementary movie 1(File)
File 2.2 MB
Supplementary material: File

De Paoli et al. supplementary movie 2

Evolution of the concentration field from simulation S4. Colour scale as shown in figure 7.
Download De Paoli et al. supplementary movie 2(File)
File 3.4 MB
Supplementary material: File

De Paoli et al. supplementary movie 3

Evolution of the concentration field from simulation S6. Colour scale as shown in figure 7.
Download De Paoli et al. supplementary movie 3(File)
File 7.5 MB
Supplementary material: File

De Paoli et al. supplementary movie 4

Evolution of the concentration field from simulation S12. Colour scale as shown in figure 7.
Download De Paoli et al. supplementary movie 4(File)
File 24.7 MB