1. Introduction
The transport and mixing of a solute or a dispersed phase within a fluid-saturated porous medium are key processes common to many subsurface flows, including petroleum migration (Simmons, Fenstemaker & Sharp Reference Simmons and Fenstemaker2001), geological sequestrations of carbon dioxide (Emami-Meybodi et al. Reference Emami-Meybodi, Hassanzadeh, Green and Ennis-King2015) and underground hydrogen storage (Krevor et al. 2023). These processes are controlled by a convective flow driven by the local fluid density differences, due to a non-uniform solute distribution within the domain. The evolution of these systems is hard to monitor due to the inaccessibility of the sites, typically located hundreds of metres underground. Nevertheless, an accurate prediction of solute transport in subsurface buoyancy-driven flows is required to address key environmental challenges, including also the design of subsurface storage sites for radioactive waste (Woods Reference Woods2015) and the development of remediation strategies for contaminated regions (LeBlanc Reference LeBlanc1984; Van Der Molen & Van Ommen Reference Van Der Molen and Van Ommen1988; Bear & Cheng Reference Bear and Cheng2010). For instance, to design groundwater remediation strategies for light and dense non-aqueous petroleum liquids, resulting from spillage of fuels or chemicals, it is key to determine the area over which these contaminants spread, and then where the remediation should take place (Woods Reference Woods2015).
Accurate predictions of the above-mentioned flows are made further complex by the interplay of convective and diffusive processes, representing the driving and dissipative mechanics, respectively, of these systems (De Paoli Reference De Paoli2023), and therefore controlling mixing. At the pore scale and in the absence of flow, solute transport is regulated by molecular diffusive mixing, acting to reduce local gradients of solute concentration: caused by the Brownian motion of the molecules, it produces a flux of solute from regions of high concentration towards regions of low concentration. When the fluid moves within the intricate interstitial pore channels, the fluid and the solute cannot penetrate into the solid obstacles (e.g. rock grains) forming the porous matrix, and will follow random-walk-type paths resulting in a further solute redistribution in the flow (longitudinal) direction (Woods Reference Woods2015). In addition, the velocity at the pore scale varies in magnitude and direction from point to point within the fluid present in the void space. This causes fluid particles travelling along their own microscopic streamline to spread out, contributing to further disperse in the flow direction any solute contained in it. The velocity gradient existing within the pores in the direction perpendicular (also defined as transverse) to the flow causes a differential transport of solute in that direction, leading to gradients of concentration, and ultimately producing solute fluxes in the transverse direction due to molecular mixing. The mechanisms of spreading produced by a flow through a porous medium are referred to as ‘mechanical’ dispersion, to indicate a spreading due to fluid mechanical phenomena. Unlike molecular diffusion which occurs also in a still fluid, both the flow and the porous medium are required for mechanical dispersion to take place. We refer the reader to Bear & Cheng (Reference Bear and Cheng2010) for an extensive presentation of this processes. We remark that molecular diffusion is a pore-scale process purely controlled by concentration gradients. Despite being originated by a very different physical mechanism, mechanical dispersion (which is a macroscopic flow property) has been modelled through a functional form analogous to that of molecular diffusion, since both contribute to spread the solute. In addition, the effects of mechanical dispersion may overcome those of molecular diffusion in the absence of flow by several orders of magnitude. At the macroscopic (Darcy) level, the combination of molecular diffusion and mechanical dispersion is also indicated as hydrodynamic dispersion (Bear & Cheng Reference Bear and Cheng2010; Wen, Chang & Hesse Reference Wen, Chang and Hesse2018). For simplicity, we refer to ‘mechanical dispersion’ as ‘dispersion’. Some of the modelling approaches proposed to describe the spreading of a solute in a porous medium include modelling of the shear (Taylor Reference Taylor1953), of the no-slip condition at the solid boundaries (Saffman Reference Saffman1959) and of the dead-end pores (Coats & Smith Reference Coats and Smith1964).
Deriving a unified model reproducing all these effects is challenging, due to the large space of the governing parameters (depending on fluid, medium and flow properties) and the different nature of the physical mechanisms involved. An approach commonly adopted consists of introducing a dispersion tensor in the transport equation governing the scalar phase (solute). This tensor combines molecular diffusion (quantified by the diffusivity coefficient
$D_m^*$
) and dispersion, modelled as an equivalent (or effective) value of solute diffusivity. Two coefficients are used to describe the spreading of solute in longitudinal (i.e. aligned with the flow) and transverse (i.e. perpendicular to the flow) directions (Delgado Reference Delgado2007), and are indicated with
$D_l^*$
and
$D_t^*$
, respectively. In general, these coefficients are obtained through laboratory experiments in constant-displacement flows through uniform porous media, but analogous results can be also achieved using pore-network models (Bijeljic, Muggeridge & Blunt Reference Bijeljic, Muggeridge and Blunt2004). It has been shown that the dispersion coefficients depend on the relative importance of flow velocity to solute diffusivity (Péclet number), among other parameters (e.g. Schmidt number and surface tension), and several flow regimes have been identified and well characterised theoretically (Saffman Reference Saffman1959; Koch & Brady Reference Koch and Brady1985; Puyguiraud, Gouze & Dentz Reference Puyguiraud, Gouze and Dentz2021).
In the above-mentioned works, the flow field is imposed a priori, and it is independent of the solute distribution within the system. In the case of buoyancy-driven flows, however, the solute field determines the flow, and the effects of dispersion have been quantified only for a few configurations. Menand & Woods (Reference Menand and Woods2005) provided measurements of the longitudinal dispersion coefficient for two miscible fluid layers with different density, and driven by either a gravitationally stable or unstable linear displacement flow, and therefore with dispersion provided by both displacement- and buoyancy-induced flows. They found that when buoyancy dominates, representing the case of interest in this work, dispersion controls the initial diffusive growth of the mixing zone and the formation and growth of the instabilities. From a macroscopic perspective, buoyancy-driven flows have been investigated by Hesse & Woods (Reference Hesse and Woods2010). They observed that the effect of dispersion associated with the presence of barriers is different from that relative to small-scale dispersion and, depending on the geometry of the barriers, may exceed the latter by several orders of magnitude. The role of dispersion in solute transport in buoyancy-driven flows has been partially explored via numerical simulations in semi-infinite systems (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; Emami-Meybodi Reference Emami-Meybodi2017; Liang et al. Reference Liang, Wen, Hesse and DiCarlo2018; Michel-Meyer et al. Reference Michel-Meyer, Shavit, Tsinober and Rosenzweig2020; Dhar et al. Reference Dhar, Meunier, Nadal and Meheust2022) and in confined steady-state flows with constant driving (Wen et al. Reference Wen, Chang and Hesse2018; Gasow et al. Reference Gasow, Kuznetsov, Avila and Jin2021). Despite these efforts, the effects of dispersion on buoyancy-driven flows in closed systems (i.e. in the absence of an external driving) remain largely unexplored, and we aim precisely at bridging this gap.
We analyse the role of dispersion in convective mixing in two-dimensional, homogeneous and isotropic porous media. We consider a Rayleigh–Taylor instability, obtained by stacking, in an unstable configuration, two fluid layers initially divided by a horizontal interface. The fluids are fully miscible, so there is no stabilising effect introduced by the presence of a surface tension, and a solute is responsible for the density differences which drive the flow. This flow configuration represents the archetypal problem used to study closed systems, and has been extensively studied in previous numerical works, in two dimensions and without dispersion (see e.g. De Wit Reference De Wit2016, Reference De Wit2020; De Paoli et al. Reference De Paoli, Zonta and Soldati2019b
; Gopalakrishnan, Knaepen & De Wit Reference Gopalakrishnan, Knaepen and De Wit2021, and references therein). Recently this problem has been also investigated in detail in three dimensions (Boffetta, Borgnino & Musacchio Reference Boffetta, Borgnino and Musacchio2020; Boffetta & Musacchio Reference Boffetta and Musacchio2022), allowing one to explore the role of the flow dimensionality. In this work, we include also the effect of dispersion of solute using the anisotropic Fickian dispersion tensor formulation proposed by Bear (Reference Bear1961). In this modelling framework, dispersion is controlled by the dimensionless parameters
$r=D_l^*/D_t^*$
, defined as the anisotropy ratio, and the relative importance of molecular diffusion to transverse dispersion,
$\varDelta =D_m^*/D_t^*$
. While
$\varDelta$
may vary over several orders of magnitude (e.g.
$0.02\leqslant \varDelta \leqslant 100$
; see Liang et al. (Reference Liang, Wen, Hesse and DiCarlo2018)),
$r$
is restricted to a narrower range, namely
$1\leqslant r\leqslant O(10)$
, where
$r\approx 10$
is appropriate for advection-dominated flows with solute transport (Bijeljic & Blunt Reference Bijeljic and Blunt2007). The relative strength of buoyancy and molecular diffusion is quantified by the Rayleigh–Darcy number (
$\textit{Ra}$
). With the aid of numerical Darcy simulations, we investigate the mixing dynamics at
$\textit{Ra}=10^4$
, in correspondence of which multiple flow regimes are present (unlike at lower values of
$\textit{Ra}$
). We employ the anisotropic dispersion tensor formulation to simulate high-
$\textit{Ra}$
convection in a closed system. We use the model proposed by Wen et al. (Reference Wen, Chang and Hesse2018), who investigated a statistically steady porous Rayleigh–Bénard system, to analyse a transient flow. Moreover, while in the Rayleigh–Bénard case heat or mass transfer at the boundaries is possible, the problem investigated here is a closed system with no external forcing, where the driving is uniquely given by the initial available potential energy. Finally, this work complements recent studies on Rayleigh–Taylor convection in porous media at high Rayleigh–Darcy numbers (De Paoli et al. Reference De Paoli, Giurgiu, Zonta and Soldati2019a
,
Reference De Paoli, Zonta and Soldatib
; Boffetta et al. Reference Boffetta, Borgnino and Musacchio2020; Boffetta & Musacchio Reference Boffetta and Musacchio2022) by systematically investigating the effect of dispersion.
The paper is organised as follows. In § 2 we present the flow configuration and the numerical method employed. The theoretical framework developed to carry out a detailed analysis of the molecular and dispersive components of the mean scalar dissipation is presented in § 3. The flow evolution is first studied in the absence of dispersion (
$\varDelta \to \infty$
; § 4), and then the effects of
$\varDelta$
(§ 5) and
$r$
(§ 6) are considered. Finally, the conclusions, an example of relevance to geophysical applications and limitations of the current approach are discussed in § 7.
2. Methodology
We study a convective flow in a porous medium at the Darcy scale. In this framework, the equations are written for quantities averaged over a reference (or representative) elementary volume (REV) (Whitaker Reference Whitaker1998), consisting of an intermediate scale larger than the individual pores, but smaller than the characteristic size of the macroscopic domain and of the flow structures. A thorough description of how to select an appropriate size of the REV is discussed by Whitaker (Reference Whitaker1998) and Nield & Bejan (Reference Nield and Bejan2006). In general, the results of the averaging procedure should be independent of the size of the REV. For buoyancy-driven flows, this procedure can be applied when: (i) the dissipative mechanisms, such as molecular diffusion and viscous dissipation, dominate over inertia (De Paoli Reference De Paoli2023); and (ii) the flow structures are large compared with the characteristic pore size (Hewitt Reference Hewitt2020).
Sketch of the flow configuration with all quantities shown in dimensionless units. An example of initial concentration field,
$C$
, consisting of a heavy fluid layer with maximum solute concentration at top (
$C=1$
) and minimum at bottom (
$C=0$
), is shown. The flow reference frame (
$x,z$
), the boundary conditions and the direction of the gravitational acceleration (
$\boldsymbol{g}$
) are indicated, as well as the domain size in horizontal (
$L$
) and vertical (
$\textit{Ra}$
) directions.
We consider a fluid-saturated porous medium in a two-dimensional domain having uniform porosity
$\phi$
and permeability
$K$
. We assume the flow is incompressible and governed by the Darcy equation, and it is characterised by an unstable density difference (
$\Delta \rho ^{*}$
) induced by the presence of a solute concentration field,
$C^*$
(the superscript
$^{*}$
is used to indicate dimensional variables). The system is periodic in the horizontal direction
$x^*$
and confined by two walls that are impermeable to fluid and solute in the vertical direction
$z^*$
, along which the gravitational acceleration
$\boldsymbol{g}$
acts. The domain extension is
$L^*_x$
and
$L^*_z$
in horizontal and vertical directions, respectively. A sketch of the computational domain with indication of the boundary conditions is reported in figure 1. The scalar field
$C^{*}$
varies between
$C^{*}_{\textit{min}}$
and
$C^{*}_{\textit{max}}$
. The evolution of this field is controlled by the advection–dispersion equation (Nield & Bejan Reference Nield and Bejan2006):
\begin{equation} \phi \frac {\partial C^{*}}{\partial t^{*}}+\boldsymbol{u}^{*}\boldsymbol{\cdot }\boldsymbol{\nabla} ^{*}C^{*}=\boldsymbol{\nabla} ^* \boldsymbol{\cdot }\left (\phi \unicode{x1D63F}^{\kern0.5pt*} \boldsymbol{\nabla} ^{*} C^{*}\right )\! \text{ ,} \end{equation}
where
$t^{*}$
is time and
$\boldsymbol{u}^{*}=(u^{*},w^{*})$
is the velocity field. The effect of dispersion is accounted for by the hydrodynamic dispersion tensor
$\unicode{x1D63F}^{\kern0.5pt*}$
that depends on the local flow conditions. We follow the formulation proposed by Bear (Reference Bear1961) and later widely employed (Bear & Bachmat Reference Bear and Bachmat2012; Emami-Meybodi Reference Emami-Meybodi2017; Wen et al. Reference Wen, Chang and Hesse2018) to model dispersion effects, and
$\unicode{x1D63F}^{\kern0.5pt*}$
reads
\begin{equation} \unicode{x1D63F}^{\kern0.5pt*} = D_m^*\unicode{x1D644}+(\alpha _l^*-\alpha _t^*)\frac {\boldsymbol{u}^*(\boldsymbol{u}^*)^{\textrm {T}}}{\lvert \boldsymbol{u}^* \rvert }+\alpha _t^*\unicode{x1D644}\lvert \boldsymbol{u}^* \rvert ,\end{equation}
where
$\unicode{x1D644}$
is the identity matrix,
$D_m^*$
the molecular diffusion coefficient,
$\alpha _l^*$
the longitudinal dispersivity and
$\alpha _t^*$
the transverse dispersivity (where
$\alpha _l^*\geqslant \alpha _t^*$
; Delgado Reference Delgado2007). Note that the advection–diffusion form of the tensor
$(\unicode{x1D63F}^{\kern0.5pt*} = D_m^*\unicode{x1D644}\,)$
is recovered when
$\boldsymbol{u}^*\alpha _l^*/D_m^*\ll 1$
.
We consider the fluid density,
$\rho ^{*}$
, to be a linear function of the concentration:
\begin{equation} \rho ^{*}(C^{*})=\rho ^{*}\left(C^{*}_{\textit{min}}\right)+\Delta \rho ^{*}\frac {C^{*}-C^{*}_{\textit{min}}}{\Delta C^*} \text{ ,} \end{equation}
with
$\Delta C^* = C^{*}_{\textit{max}}-C^{*}_{\textit{min}}$
and
$\Delta \rho ^{*}=\rho ^{*}(C^{*}_{\textit{max}})-\rho ^{*}(C^{*}_{\textit{min}})$
. Assuming the validity of the Boussinesq approximation, which is reasonable, e.g. in the context of brine transport in porous media (Landman & Schotting Reference Landman and Schotting2007), the flow field is fully described by the continuity and the Darcy equations:
\begin{align} \boldsymbol{\nabla} ^{*}\boldsymbol{\cdot }\boldsymbol{u}^{*}=0,\quad \boldsymbol{u}^{*}=-\frac {K}{\mu }\left (\boldsymbol{\nabla} ^{*} P^{*}+\rho ^{*}g \boldsymbol{k}\right )\!,\end{align}
with
$\mu$
the fluid viscosity (constant),
$P^{*}$
the pressure and
$\boldsymbol{k}$
the vertical unit vector. Additional correction terms to account for inertial effects can be included (Nield & Bejan Reference Nield and Bejan2006), but lie beyond the scope of this work. Since the walls are impermeable to the fluid, the boundary condition reads
\begin{equation} \boldsymbol{u}^* \boldsymbol{\cdot }\boldsymbol{n} = 0 \, \Rightarrow \, \begin{cases} w^*\left(z^*=-L^{*}_{z}/2\right)=0\\ w^*\left(z^*=+L^{*}_{z}/2\right)=0 \end{cases} \end{equation}
with
$\boldsymbol{n}$
the unit vector perpendicular to the boundary (note that slip at the walls is possible). At the upper and lower walls, no-flux (
$\partial _{z^*} C^*=0$
) conditions are considered. Periodicity is forced in the wall-parallel directions.
2.1. Dimensionless equations
A natural velocity scale relevant to the convective system examined is the buoyancy velocity,
$\mathcal{U}^{*}=g \Delta \rho ^{*} K / \mu$
. After the onset of convection, fingers start to develop vertically from the centreline (
$z^*=0$
), and the domain height is relevant only after the fingers reach the walls. Therefore, it may be convenient to make the equations dimensionless with respect to flow units that are independent of the domain geometry. In particular, as proposed by Fu, Cueto-Felgueroso & Juanes (Reference Fu, Cueto-Felgueroso and Juanes2013), one can use as a reference length scale
$\ell ^*=\phi D_m^*/\mathcal{U}^*$
, where the vertical domain extension
$L^*_z$
is not included. Using these scales we obtain the following set of dimensionless variables:
\begin{align} C&=\frac {C^{*}-C^{*}_{\textit{min}}}{C^{*}_{\textit{max}}-C^{*}_{\textit{min}}},\quad x=\frac {x^{*}}{\ell ^{*}},\quad \boldsymbol{u}=\frac {\boldsymbol{u}^{*}}{\mathcal{U}^{*}}, \end{align}
\begin{align} t&=\frac {t^{*}}{\phi \ell ^{*}/\mathcal{U}^{*}},\quad p=\frac {p^{*}}{\Delta \rho ^{*}g\ell ^{*}} , \end{align}
where we introduced the reduced pressure
$p^{*}=P^*+\rho ^*(C^*_{\textit{min}})gz^*$
. We finally derive the dimensionless form of the governing equations (2.1)–(2.4):
\begin{align}& \frac {\partial C}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }C = \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \unicode{x1D63F} \boldsymbol{\nabla }C\right )\!, \end{align}
\begin{align}&\qquad\quad\,\, \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{align}
\begin{align}&\qquad \boldsymbol{u}=-\left ( \boldsymbol{\nabla }p + C \boldsymbol{k}\right )\! , \end{align}
where
\begin{equation} \textit{Ra}=\frac {g \Delta \rho ^{*} K L^{*}_{z} }{ \phi D_m^* \mu }=\frac {\mathcal{U}^{*} L^{*}_{z}}{\phi D_m^*} \end{equation}
is the Rayleigh–Darcy number (indicated as the Rayleigh number in the following). We recall here that
$K$
and
$\phi$
represent the medium permeability and porosity, respectively, and are considered uniform here. Finally, the dispersion tensor
$\unicode{x1D63F}^{\kern0.5pt*}$
introduced in (2.2) and expressed in dimensionless form reads
\begin{equation} \unicode{x1D63F}=\unicode{x1D644}+\frac {1}{\varDelta }\left [(r-1)\frac {\boldsymbol{u}\boldsymbol{u}^{\textrm {T}}}{\lvert \boldsymbol{u}\rvert }+\unicode{x1D644}\lvert \boldsymbol{u}\rvert \right ]\! , \end{equation}
where
\begin{equation} \varDelta =\frac {D_m^*}{D_t^*},\quad r=\frac {D_l^*}{D_t^*}=\frac {\alpha _l^*}{\alpha _t^*} \end{equation}
with
$D_m^*$
the molecular diffusion coefficient,
$D_t^*=\alpha _t^*\mathcal{U}^*$
the transverse dispersion coefficient,
$D_l^*=\alpha _l^*\mathcal{U}^*$
the longitudinal dispersion coefficient and
$r$
the dispersivity ratio. The flow is completely defined by four dimensionless parameters:
$\textit{Ra}$
,
$L=L^*_x/\ell ^*$
,
$\varDelta$
and
$r$
. This formulation is similar to that proposed by Wen et al. (Reference Wen, Chang and Hesse2018). However, in this case the governing parameter
$\textit{Ra}$
does not appear explicitly in the equations, and it corresponds to the dimensionless domain height (
$\textit{Ra} = L^{*}_{z} /\ell ^*$
; see also figure 1). In addition, provided
$L$
is sufficiently large to allow the formation of multiple flow structures and to minimise the effects of the periodic forcing, the flow is also independent of
$L$
itself.
2.2. Numerical solution of the equations
The set of equations (2.8)–(2.10) is solved numerically with the aid of the second-order finite-difference code AFiD-Darcy open-sourced by our research group (De Paoli Reference De Paoli2025). The code consists of an efficient solver for massively parallel simulations of convective, wall-bounded and incompressible porous media flows, and it is based on the initial version of AFiD developed for turbulent flows (Van Der Poel et al. Reference Van Der Poel, Ostilla-Mónico, Donners and Verzicco2015). The algorithm is based on a pressure-correction scheme and employs an efficient fast-Fourier-transform-based solver. The parallelisation method is implemented in a two-dimensional pencil-like domain decomposition, which enables efficient parallel large-scale simulations. The implementation in the absence of dispersion was validated against different canonical flows in porous media, including Rayleigh–Taylor instability (Boffetta et al. Reference Boffetta, Borgnino and Musacchio2020). Additional numerical details of the solution algorithm are provided by De Paoli et al. (2025a). In addition, in this work we solve the flow including the effect of dispersion (2.12) (see Appendix B.1 for further details on the numerical treatment of the dispersive terms). The algorithm with dispersion is validated against the work of Wen et al. (Reference Wen, Chang and Hesse2018) repeating their simulations in Rayleigh–Bénard–Darcy configuration, and comparing the results in terms of molecular, dispersive and total Nusselt numbers.
We employ here a uniform grid spacing with cells of size
$\Delta x\approx \Delta z\lt 10$
(expressed in dimensionless units as discussed in § 2.1) in horizontal and vertical directions, respectively. The present simulations are over-resolved compared with the minimum requirements of the case without dispersion (corresponding to
$\Delta x = \Delta z = 15.625$
; see De Paoli et al. Reference De Paoli2025a
). The minimum number of grid points in the vertical direction used here is 128, making the grid even further over-resolved in the low-
$\textit{Ra}$
case. In the horizontal direction, the domain size is fixed (
$L=10^5$
), as well as its resolution (
$N_x=10\,240$
,
$\Delta x = L/N_x\lt 10$
), for all simulations except for the smallest
$\varDelta$
considered here (
$\varDelta =5\times 10^{-2}$
), for which the resolution is doubled and the domain width halved in order to keep the computational cost affordable. The details of all simulations performed are listed in table 1. The flow is initialised with a step-like profile obtained from the analytical diffusive solution (2.8). Additional details of the initial condition and the grid requirements are also provided in Appendix B.2.
Summary of the parameters employed in the simulations. The governing parameters of the flow (Rayleigh number
$\textit{Ra}$
, domain width
$L$
and domain aspect ratio
$L/\textit{Ra}$
) and the dispersion parameters (
$\varDelta$
and
$r$
; see (2.13)) are indicated, as well as the grid resolution employed. Finally, the dimensionless growth rate of the mixing layer,
$\gamma$
, defined in (4.4), is reported for the simulations with
$\textit{Ra}=10^4$
.
3. Global budgets and mixing indicators
Exact global conservation equations relative to the transport of the scalar quantity can be derived and employed to investigate the mixing process. We consider here the case of an incompressible and dispersive flow governed by (2.8)–(2.12), and we derive the first- and second-order global budgets in §§ 3.1 and 3.2, respectively. Finally, we introduce the degree of mixing, representative of the current mixing state of the system, in § 3.3.
3.1. First-order global budget
We apply the volume-average operator
$\langle \boldsymbol{\cdot }\rangle =1/V\int _V\boldsymbol{\cdot }\textrm { d}V$
to (2.8). Using the divergence theorem, the incompressibility of the flow and the no-penetration/periodic boundary conditions, we derive the global flux of solute in the domain:
\begin{align} \frac {{\rm d} \langle C\rangle }{{\rm d} t} &= \frac {1}{V}\int _{S=\partial V} \left [ \left (\unicode{x1D63F}\boldsymbol{\nabla }C \right )_{z=-\textit{Ra}/2} +\left (\unicode{x1D63F}\boldsymbol{\nabla }C \right )_{z=+\textit{Ra}/2} \right ]\boldsymbol{\cdot }\boldsymbol{n}\textrm { d}S \end{align}
\begin{align} &= \frac {1}{\textit{Ra} S}\int _{S} \left [ -\left (D_{zz}\frac {\partial C}{\partial z} \right )_{z=-\textit{Ra}/2} +\left (D_{zz}\frac {\partial C}{\partial z} \right )_{z=+\textit{Ra}/2} \right ]\textrm { d}S \end{align}
\begin{align} &= \frac {1}{\textit{Ra}}\Bigl [F\left (z=\textit{Ra}/2\right )-F\left (z=-\textit{Ra}/2\right )\Bigr ] , \end{align}
where
$D_{zz} = 1+\varDelta ^{-1} [\lvert \boldsymbol{u}\rvert +(r-1)w^{2}/\lvert \boldsymbol{u}\rvert ]$
(see Appendix A) and
$F(z_i)$
is the dispersive flux at the walls, defined as
\begin{equation} F(z_i) = \frac {1}{S}\int _S\left [\left (\frac {\partial C}{\partial z}\right )_{z=z_i} +\left (\frac {|\boldsymbol{u}|}{\Delta }\frac {\partial C}{\partial z}\right )_{z=z_i}\right ]\textrm {d}S \end{equation}
with
$S=L_x$
or
$S=L_x\times L_y$
in two and three dimensions, respectively. Note that (3.4) corresponds to the definition of (10) proposed by Wen et al. (Reference Wen, Chang and Hesse2018). In the present configuration, due to the no-flux condition at the boundaries (
$\partial _z C=0$
), we have that
$F=0$
and then
\begin{equation} \frac {{\rm d} \langle C\rangle }{{\rm d}t} = 0 . \end{equation}
3.2. Second-order global budget
Following previous works on convection in semi-infinite layers (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012), Rayleigh–Bénard flows (Otero et al. Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004; Hassanzadeh; Chini & Doering Reference Hassanzadeh, Chini and Doering2014; Zhu, Fu & De Paoli Reference Zhu, Fu and De Paoli2024; Hu & Yang Reference Hu and Yang2024) and Hele-Shaw flows (Letelier, Mujica & Ortega Reference Letelier, Mujica and Ortega2019; Ulloa & Letelier Reference Ulloa and Letelier2022; Ulloa, Noto & Letelier Reference Ulloa, Noto and Letelier2025), we multiply (2.8) by
$C$
and apply the volume-average operator
$\langle \boldsymbol{\cdot }\rangle$
. We use the divergence theorem and the incompressibility of the flow to derive the global budget (De Paoli Reference De Paoli2023):
\begin{align} \frac {1}{2}\frac {{\rm d} \langle C^2\rangle }{{\rm d}t} &= \frac {1}{V}\int _S \left [ \left (C\unicode{x1D63F}\boldsymbol{\nabla }C \right )_{z=-\textit{Ra}/2} +\left (C\unicode{x1D63F}\boldsymbol{\nabla }C \right )_{z=+\textit{Ra}/2} \right ]\boldsymbol{\cdot }\boldsymbol{n}\textrm { d}S- \langle (\boldsymbol{\nabla }C)\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C)\rangle \end{align}
\begin{align} &=\frac {1}{\textit{Ra}}\Bigl [-(CF)_{z=-\textit{Ra}/2}+(CF)_{z=+\textit{Ra}/2}\Bigr ]-\langle (\boldsymbol{\nabla }C)\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C)\rangle \end{align}
with
$\boldsymbol{n}$
a unit vector normal to the boundary. Since we refer to a Rayleigh–Taylor configuration, i.e.
$\partial _z C(z=\pm \textit{Ra}/2)=0$
, we have that
\begin{equation} \frac {1}{2}\frac {{\rm d} \langle C^2\rangle }{{\rm d}t} = -\langle (\boldsymbol{\nabla }C)\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C)\rangle . \end{equation}
Inspired by the Grossmann–Lohse theory (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Lohse & Shishkina Reference Lohse and Shishkina2024) where the key idea is to spatially split the viscous dissipation rate into boundary-layer and bulk contributions, here we divide the scalar dissipation into two parts, one ascribed to the molecular component and the other to the dispersive component of the dispersion tensor. Expressed in terms of mean scalar dissipation, (3.8) reads
\begin{equation} \frac {1}{2}\frac {{\rm d} \langle C^2\rangle }{{\rm d}t} =-\frac {1}{\textit{Ra}}\chi , \end{equation}
where
\begin{align} \chi = \chi _m + \chi _d ,\quad \chi _m=\textit{Ra}\langle |\boldsymbol{\nabla }C|^2\rangle ,\quad \chi _d=\textit{Ra}\big\langle (\boldsymbol{\nabla }C)\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C)-|\boldsymbol{\nabla }C|^2\big\rangle. \end{align}
The components of the mean scalar dissipation introduced in (3.10) are defined as the total, molecular and dispersive contributions, respectively. Note that in the absence of dispersion (
$\varDelta \to \infty$
) we obtain
$\unicode{x1D63F}=D_m^*\unicode{x1D644}$
and then
$\chi _d=0$
. The choice to introduce
$\textit{Ra}$
in the definition of mean dissipation is required to obtain results that are self-similar when different
$\textit{Ra}$
are considered, i.e. to make the results independent of the domain height until the fingers reach the horizontal boundaries (De Paoli et al. Reference De Paoli, Giurgiu, Zonta and Soldati2019a
). In dimensional terms, the components of the mean dissipation defined in (3.10) read
\begin{align} \chi ^* &= \chi ^*_m + \chi ^*_d ,\quad \chi ^*_m = D_m^*\langle |\boldsymbol{\nabla} ^* C^*|^2\rangle =\frac {\mathcal{U}^*(\Delta C^*)^2}{\phi H^*}\chi _m ,\end{align}
\begin{align} \chi ^*_d &= \langle (\boldsymbol{\nabla} ^* C^*)\boldsymbol{\cdot }(\unicode{x1D63F}^{\kern0.5pt*}\boldsymbol{\nabla} ^* C^*)-|\boldsymbol{\nabla} ^* C^*|^2\rangle = \frac {\mathcal{U}^*(\Delta C^*)^2}{\phi H^*}\chi _d. \end{align}
Note that
$\chi _m^*$
in (3.11) matches the definition of dissipation used by De Paoli et al. (Reference De Paoli, Howland, Verzicco and Lohse2024).
3.3. Quantification of mixing
Following the approach proposed by Jha, Cueto-Felgueroso & Juanes (Reference Jha, Cueto-Felgueroso and Juanes2011), we introduce the degree of mixing, a quantity that varies between 0 and 1, and is representative of the current mixing state of the flow. The initial state of the system is characterised by two uniform layers having different concentration and divided by a sharp interface, i.e.
$C(x,z\gt 0,t=0)=1$
and
$C(x,z\lt 0,t=0)=0$
. Ultimately, the system will achieve a well-mixed condition and the concentration field will be uniform, corresponding to
$C(x,z,t\to \infty )=1/2$
. As a result, the mean concentration variance, defined as
$\sigma ^2=\langle C^2 \rangle - \langle C \rangle ^2$
, is initially maximum and equal to
$\sigma ^2(t=0)=\sigma ^2_{\textit{max}}=1/4$
, and ultimately minimum and corresponding to
$\sigma ^2(t\to \infty )=0$
. The variance is representative of the mixing state, which we quantify using the degree of mixing:
\begin{equation} M(t) = 1-\frac {\sigma ^2(t)}{\sigma ^2_{\textit{max}}}. \end{equation}
Such definition gives
$M(t=0)=0$
when the two layers are perfectly segregated and
$M(t\to \infty )=1$
when a perfect mixing is achieved. Also in this case, we split the degree of mixing into two contributions resulting from the partition of the dissipation (Grossmann & Lohse Reference Grossmann and Lohse2000,Reference Grossmann and Lohse2001). Given the definition of
$\sigma ^2$
and the budget (3.9), one can rewrite (3.13) as
\begin{align} M(t) &= 1-\frac {\sigma ^2(t)}{\sigma ^2_{\textit{max}}}\nonumber \\ &=1-\frac { \langle C^2(0) \rangle -\langle C(t) \rangle ^2 + \int _{\langle C^2(0)\rangle }^{\langle C^2(t)\rangle } \text{ d} \langle C^2\rangle }{\sigma ^2_{\textit{max}}}\nonumber \\ &=\frac {2}{\sigma ^2_{\textit{max}}\textit{Ra}}\int _0^t\left (\chi _m + \chi _d\right ) \text{ d}\tau = M_m(t) + M_d(t), \end{align}
where the degree of mixing has been split into the molecular and dispersive components, respectively:
\begin{align} M_m(t) = \frac {2}{\sigma ^2_{\textit{max}}\textit{Ra}}\int _0^t \chi _m \text{ d}\tau ,\quad M_d(t)=\frac {2}{\sigma ^2_{\textit{max}}\textit{Ra}}\int _0^t \chi _d \text{ d}\tau . \end{align}
4. Flow evolution without dispersion (
$\boldsymbol{\varDelta \to \infty}$
)
We consider the simulations performed without dispersion corresponding to
$\varDelta \to \infty$
in (2.12), and we explore the dynamics for different values of Rayleigh number
$\textit{Ra}$
. We cover a wide range of
$\textit{Ra}$
, namely
$10^2\leqslant \textit{Ra} \leqslant 2\times 10^4$
, and we simulate domains having constant width
$L=10^5$
. The details of all simulations performed are listed in table 1. In the absence of dispersion, the influence of
$\textit{Ra}$
on the flow dynamics has been previously investigated in several two-dimensional works (De Wit Reference De Wit2004; De Paoli et al. Reference De Paoli, Giurgiu, Zonta and Soldati2019a
,
Reference De Paoli, Zonta and Soldatib
; Borgnino, Boffetta & Musacchio Reference Borgnino, Boffetta and Musacchio2021). Here we complement these analyses, which are essential to derive a clear picture of the flow dynamics also in the case of dispersive flows. For all
$\textit{Ra}$
considered, we report in figure 2 the evolution of the molecular component of the mean scalar dissipation
$\chi _m$
(for
$\varDelta \to \infty$
,
$\chi _d=0$
). The dimensionless set of variables used is particularly suitable to highlight the self-similar behaviour of the system. In particular, the flow evolution is independent of
$\textit{Ra}$
until the domain walls have an effect on the flow, i.e. until the fingers grow and their flow field feels the influence of the horizontal impermeable walls. For very small values of
$\textit{Ra}$
(
$\textit{Ra}\leqslant 10^3$
), the dynamics is soon influenced by the presence of the walls, and convective instabilities may not form or only partially develop.
Evolution of the molecular mean scalar dissipation
$\chi _m$
in the absence of dispersion and for different Rayleigh numbers
$\textit{Ra}$
(solid lines). The flow evolution is independent of
$\textit{Ra}$
until the flow field of the fingers is significantly influenced by the presence of the horizontal walls. The analytical diffusive solutions in the unconfined (4.2) (dotted line) and confined (C5) (dashed line) cases are also reported.
Evolution of the concentration field relative to the simulation
$\textit{Ra}=10^4$
. (a–e) A portion of the domain is shown, corresponding to half of the domain width (left panels, indicated with (i)) and 1/20 of the domain width (right panels, indicated with (ii) and corresponding to the white rectangle in the corresponding panels (i)). The entire domain simulated is shown in (d–h). The data correspond to the points indicated in figure 4.
Evolution of the mean scalar dissipation for the simulation
$\textit{Ra}=10^4$
without dispersion (see table 1 for additional details). The concentration fields and profiles correspond to the instants indicated by the letters shown in figure 3. The diffusive solution (4.2) is also reported in (a). The flow regimes identified are indicated in (b).
In the following, to analyse in detail the flow dynamics, we focus on the simulation at
$\textit{Ra}=10^4$
, which serves also as a reference for the dispersive cases studied in §§ 5 and 6. The system behaviour relative to different
$\textit{Ra}$
can be inferred from the analysis presented. The flow evolution is initially characterised by inspecting the concentration fields and the concentration profiles in figure 3, together with the corresponding mixing indicators in figure 4. The process is also illustrated in supplementary movie S1 available at https://doi.org/10.1017/jfm.2025.10684. We identified several flow regimes characterising the evolution of the system, which we discuss individually in the following.
The flow is initialised with a step-like concentration profile (figure 3
a) with a thin interface dividing the two fluid layers and no flows (
$\boldsymbol{u}=0$
). Therefore, the interface will initially grow driven by diffusion, making the concentration gradient (and then the mean scalar dissipation
$\chi =\chi _m$
; see figure 4
a) progressively reduce. This trend continues until the first convective instabilities appear at
$t\approx 10^3$
; see figure 3(b-iii). Similarly to what has been previously observed in convection in semi-infinite domains (Elenius & Johannsen Reference Elenius and Johannsen2012), the precise time at which this growth starts, as well as the maximum values of scalar dissipation later achieved, depend on the initial perturbation. However, the overall dynamics is qualitatively comparable to that described here for the initial condition considered discussed in detail in Appendix B.2. The flow is initially purely controlled by diffusion. An analytical solution of the advection–dispersion equation (2.8) can be obtained for an unconfined domain, in the absence of convection (
$\boldsymbol{u}=0$
) and assuming homogeneity in wall-parallel directions (
$\partial _x \boldsymbol{\cdot }= 0$
), which reads (De Paoli et al. Reference De Paoli, Zonta and Soldati2019b
)
\begin{equation} C(z,t)=\frac {1}{2}\left [1+\text{erf}\left (\frac {z}{2\sqrt {t}}\right )\right ]\!. \end{equation}
Results relative to the horizontally averaged concentration profiles reported in figure 5(b) indicate an excellent agreement with (4.1) (dotted line). Employing (4.1) in the definition (3.10), it follows that the contributions of the mean scalar dissipation during the initial phase (assuming
$\boldsymbol{u}=0$
) are
\begin{align} \chi _m = \frac {1}{\sqrt {8\pi t}}, \quad \chi _d = 0, \quad \chi = \frac {1}{\sqrt {8\pi t}}, \end{align}
and they are very well captured by the simulation (see figure 4 a). Correspondingly, (3.15) predicts the degree of mixing to evolve as
\begin{align} M_m =\frac {8\sqrt {t}}{\sqrt {2\pi }} , \quad M_d = 0 , \quad M = \frac {8\sqrt {t}}{\sqrt {2\pi }}. \end{align}
Note that for very low Rayleigh numbers (
$\textit{Ra}=10^2$
), the solution (4.2) does not represent a good approximation due to the effect of confinement: a no-flux boundary condition should be employed to determine the analytical solution. In Appendix C.1, we derive the analytical diffusive solution (C5) relative to the confined case, which we report in figure 2 (dashed line), which describes very well the evolution observed numerically.
Similarly to what has been observed in semi-infinite domains (see Slim (Reference Slim2014); De Paoli et al. (Reference De Paoli2025b
) and references therein), the flow is later characterised by a dissipation growth (figures 3
b,c). When the newly formed fingers have grown, they accelerate vertically sharpening the concentration gradient at their interface, corresponding to the increase in the mean dissipation observed in figure 4(a) between
$t\approx 10^3$
and
$t\approx 2.5\times 10^3$
.
Later, fingers grow further, but in an uneven manner (figures 3 c,d): some fingers extend vertically at large velocity and invade regions of the domain with uniform concentration, and far from other fingers. In this situation, the strong concentration difference between the inner and outer finger region, combined with the absence of neighbouring flow structures (and then no strong shear) prevents the fingers from expanding laterally, and makes the interface to grow diffusively while the fingers continue to extend vertically. As a result, the concentration gradients reduce and so does and the mean scalar dissipation in figure 4(a).
The velocity field generated by each finger induces interactions with the neighbouring ones, eventually leading to merging: The pioneering fingers perturb the velocity field of the neighbouring shorter fingers, compressing and forcing them to retreat and eventually to merge (figure 3 d,e). In this phase, similarly to what has been discussed in the uneven finger growth regime, the fingers keep growing vertically in unexplored regions of the flow, but now there is a large and nearly uniform reservoir of solute at the core of the fingers. The gradient of concentration at the finger interface still decreases, but at a lower rate compared with the increase of interfacial extension of the fingers, driven by their vertical growth. As a result, the overall behaviour corresponds to an increase of the mean scalar dissipation in figure 4(a).
Evolution of the horizontally averaged concentration profiles,
$\overline {C}$
, relative to simulation
$\textit{Ra}=10^4$
,
$\varDelta \to \infty$
. Profiles reported correspond to instants taken preceding the finger impact on the walls. Specifically, they are in the range
$400\leqslant t\leqslant 16\,000$
(a), in the initial diffusive regime (
$t\leqslant 1.4\times 10^3$
) (b) and in the finger merging and growth regime (
$t\geqslant 7\times 10^3$
) (c). In (b), the dashed line indicates the initial diffusive solution (4.1). In (c), the wall-normal coordinate is rescaled with
$t-t_0$
, where
$t_0=4\times 10^3$
. The dashed line represents (4.4).
During the finger merging and growth phase, the mixing region grows approximately linearly in time (namely,
$\sim\!t^{1.2}$
; see De Paoli et al. Reference De Paoli, Zonta and Soldati2019b
, Reference De Paoli, Perissutti, Marchioli and Soldati2022). Indeed, rescaling the profiles in figure 5(c) with respect to
$z/t$
, they result in being well approximated by the linear function
\begin{equation} \overline {C} = \frac {1}{2}+\frac {z}{\gamma (t-t_0)} \end{equation}
(dashed line in figure 5
c), where
$t_0=4\times 10^3$
corresponds (approximately) to the time at which this regimes starts and
$\gamma =0.59$
is the dimensionless growth rate of the mixing layer, in good agreement with the measurements of Boffetta & Musacchio (Reference Boffetta and Musacchio2022), who reported a value of 0.67 obtained at larger
$\textit{Ra}$
. The value of
$\gamma$
computed for all the simulations is reported in table 1. It has been obtained as a best fit of (4.4) for the concentration profile within the range
$0.05\leqslant C\leqslant 0.95$
, with
$t_0=4\times 10^3$
and within the interval
$7\times 10^3 \leqslant t \leqslant 1.6\times 10^4$
. For simulations
$\varDelta =5\times 10^{-2}$
and
$r=20$
, in which the finger development is slower, the profiles are considered in the interval
$1.5\times 10^4 \leqslant t \leqslant 2.7\times 10^4$
.
When the fingers reach the horizontal boundaries of the domain (
$t\approx 2\times 10^4$
), the growth of mean scalar dissipation arrests. At this time, the concentration profile is nearly linear and still unstable (figure 3
e-ii), indicating that the upper portion of the system is still characterised by a higher concentration of solute, and then larger density, compared with the lower layer. The density contrast between the upper and lower domains, however, is much smaller than at the beginning. There is no further fresh fluid available to mix with the fluid at the core of the fingers, and the local concentration gradient across the interface of the fingers decreases progressively. At the same time, the situation at the domain centre (
$z=0$
) is steady, with no new fingers forming and no fingers merging, and with the solute being transported along the finger path already created. This process occurs from the unstably stratified configuration in figure 3(e) until a stably stratified solute distribution is achieved (figure 3
g). Note that the footprint of the fingers persists for a long time (figure 3
f), and it smoothly disappears when the residual convective driving dies out. Indeed, despite being separated by approximately
$1.2\times 10^5$
dimensionless time units, the difference between the horizontally averaged concentration profiles in figures 3( f) and 3(g) is very limited.
When the stably stratified distribution (figure 3
g-ii) is achieved, the buoyancy effects are negligible, and diffusion drives mixing from now on. The system enters a final diffusive regime that, despite being characterised by a very low value of the mean dissipation (figure 4
a), is responsible for approximately 14 % of the overall mixing, i.e. (see figure 4
b)
$M$
grows from
$\approx 0.86$
in (g) to
$\approx 1$
in (h). Similarly to what we have derived for
$\textit{Ra}=10^2$
for the initial confined diffusive evolution (see Appendix C.1), we can derive here an analytical solution for the mean scalar dissipation during the late stage. Again, we use (2.8), assume that the fluid is still (
$\boldsymbol{u}=0$
) and consider the system uniform in the horizontal direction (
$\partial _x C = 0$
). At large
$\textit{Ra}$
a different initial condition has to be considered, namely a linear concentration profile, representative of a stably stratified flow that starts at
$t=t_f$
(see Appendix C.2 for addition details). The solution (C7) computed using
$n=100$
is shown in figure 6(a) for three values of
$\textit{Ra}$
, namely
$\textit{Ra}=10^2$
,
$10^3$
and
$10^4$
(with the case
$10^4$
previously discussed throughout the entire flow evolution). The accuracy of the analytical solution in predicting the decay of dissipation obtained in the simulations is excellent. As already observed in the initial confined diffusion for
$\textit{Ra}=10^2$
, also in (C7) we have that
$\partial C/\partial \bar {z}\sim \exp (-n^2)$
, indicating a fast decay with
$n$
. Therefore, to leading order the dissipation can be approximated as
\begin{equation} \chi _m(t)\sim \exp {\left [-2\left (\frac {\pi }{\textit{Ra}}\right )^2t\right ]}, \end{equation}
which we show in figure 6(a) and indicate with
$n=1$
. Also in this case, the agreement is excellent, and we conclude that (4.5) can be employed to approximate the behaviour during the final diffusive phase. For long times, the well-mixed condition is achieved, corresponding to the asymptotic equilibrium profile:
\begin{equation} C(\bar {z},t\to \infty ) = 1/2, \end{equation}
with the mean scalar dissipation (C7) and degree of mixing given by
\begin{align} \chi _m(t\to \infty )=0\quad \text{ and}\quad M_m(t\to \infty )=1, \end{align}
respectively. One can indeed observe in figure 6(b) that the well-mixed condition is ultimately achieved. Note that, for ease of comparison, the degree of mixing
$M_m$
is shown rescaled by
$\textit{Ra}$
, and the corresponding diffusive solution is
$M_m\times \textit{Ra} = 8\textit{Ra}\sqrt {t/(2\pi )}$
.
Scalar dissipation rate (a) and rescaled degree of mixing (b) are reported as a function of time
$t$
for three values of Rayleigh number,
$\textit{Ra}$
, namely
$10^2,\ 10^3,\ 10^4$
. The system is self-similar and at early times it follows the analytical solutions (initial unconfined diffusion), (4.2) for
$\chi _m$
and (4.3) for
$M_m$
, indicated here with black dotted lines. As soon as the system achieves a stably stratified condition, it enters the final diffusive phase. The scalar dissipation evolves according to (C7) indicated in (a) by the black solid lines and computed using
$n=100$
(note that for
$\textit{Ra}=10^2$
, the initial confined solution (C5) is used). However, these solutions are very well approximated also when
$n=1$
, scales as (4.5) and is indicated by the blue dashed lines. Ultimately, the domain attains the fully mixed condition (4.7) (black dashed line in b).
5. Flow evolution with dispersion: influence of
$\boldsymbol{\varDelta}$
We present here the results relating to the effect of
$\varDelta$
at
$\textit{Ra}=10^4$
and
$r=10$
. The choice of
$\textit{Ra}$
is motivated by the fact that at such
$\textit{Ra}$
all the phases of the flow evolution described in § 4 can be observed. Since the dispersion model (2.12) considered in this work relies on two parameters,
$\varDelta$
and
$r$
defined in (2.13), to keep the computational costs affordable, we decided to vary
$\varDelta$
and fix the anisotropy ratio to
$r=10$
. Determining which value is more appropriate for
$r$
is a matter of debate (Delgado Reference Delgado2007), since the longitudinal (
$D_l^*$
) and transverse (
$D_t^*$
) components of dispersion depend on several flow parameters (e.g. Schmidt number, Reynolds number, tortuosity of the medium, Péclet number, fluid phases involved; see De Paoli (Reference De Paoli2023) and references therein). The value of
$r$
varies with the dispersion regime, and in general it is not appropriate to take
$D_t^*$
to be one order of magnitude smaller than
$D_l^*$
, which is the most common choice in numerical studies of dispersion. However, when studying solute transport in the advection-dominated regime, this assumption holds (Bijeljic & Blunt Reference Bijeljic and Blunt2007), and therefore we consider our findings representative of this case.
5.1. Flow dynamics
In the presence of dispersion, the evolution of the flow follows a behaviour that is similar to the case discussed in § 4, and eventually with fewer regimes as
$\varDelta$
is decreased (see supplementary movie S2 for the time-dependent evolution of the simulation with
$\varDelta =10^{-1}$
, reported in figure 7
e). The concentration distribution over the entire field at
$t=2\times 10^4$
is reported in figure 7 for different values of
$\varDelta$
. The effects of
$\varDelta$
on the flow morphology, which are discussed in the following, are multiple.
Concentration fields at
$t=2\times 10^4$
for (a–f) different values of the dispersion parameter
$\varDelta$
. The field (a) corresponds to the case without dispersion (
$\varDelta \to \infty$
). See supplementary movie S2 for the time-dependent evolution of the simulation with
$\varDelta =10^{-1}$
(e).
Evolution of the horizontally averaged concentration profiles,
$\overline {C}$
, relative to simulation
$\textit{Ra}=10^4$
,
$\varDelta =0.1$
and
$r=10$
. Profiles reported correspond to instants taken preceding the finger impact on the walls. Specifically, they are in the range
$400\leqslant t\leqslant 16000$
(a), in the initial diffusive regime (
$t\leqslant 1.4\times 10^3$
) (b) and in the finger merging and growth regime (
$t\geqslant 7\times 10^3$
) (c). In (b), the dashed line indicates the initial diffusive solution (4.1). In (c), the wall-normal coordinate is rescaled with
$t-t_0$
, where
$t_0=4\times 10^3$
. The dashed line represents (4.4).
Results relating to the horizontally averaged concentration profiles are reported in figure 8, where different times, namely
$400\leqslant t\leqslant 16\,000$
, are shown for the simulation with
$\varDelta =0.1$
and
$r=10$
(figure 8
a). At early times (figure 8
b), the profiles follow the diffusive growth ((4.1), dashed line). Later, the mixing region extends approximately linearly in time and the profiles are very well fitted by the linear function (4.4) (dashed line in figure 8
c), where
$t_0=4\times 10^3$
and
$\gamma =0.49$
. Note that in this case the dimensionless growth of the mixing layer
$\gamma$
is lower compared with the case without dispersion (
$\gamma =0.59$
). This can be easily verified comparing the concentration fields of figures 7(a) and 7(e).
The distributions of (a) molecular (
$\textit{Ra}|\boldsymbol{\nabla }C|^2$
, see (3.10)) and (b) dispersive (
$\textit{Ra}[ (\boldsymbol{\nabla }C)\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C)-|\boldsymbol{\nabla }C|^2]$
) scalar dissipation taken at time
$t=10^4$
. The case without dispersion is shown in (i) (
$\textit{Ra}=10^4,\varDelta \to \infty$
) and the case with dispersion in (ii) (
$\textit{Ra}=10^4,\ \varDelta =10^{-1},\ r=10$
). For better visualisation, only a small region in the core of the domain is shown (
$0\leqslant x/\textit{Ra} \leqslant 1/2, -1/4\leqslant z/\textit{Ra} \leqslant 1/4$
). The dashed lines represent the iso-contours
$C=1/4$
and
$C=3/4$
. (c) The probability density function (p.d.f.) of the components of the molecular, dispersive and total dissipation normalised by their respective root mean squares (r.m.s.) and relative to the specific fields considered. For better visualisation the limits of the colourbars in (a,b) are reduced compared with the maximum/minimum values present in the field.
A first macroscopic observation is that reducing
$\varDelta$
has the effect of diminishing the number of fingers. This is apparent for
$\varDelta \leqslant 1$
(figure 7
d–f), but it occurs also for
$\varDelta =10$
(figure 7
c), while no difference is observed when
$\varDelta =10^5$
(figure 7
b) compared with the case without dispersion (
$\varDelta \to \infty$
; figure 7
a). Since for small
$\varDelta$
the neighbouring fingers are further away compared with the case without dispersion, each finger has more space available to diffusively grow horizontally. This morphological effect combines with the additional spreading provided by dispersion, which further reduces the concentration gradient across the finger interface. A visual interpretation of this mechanisms is reported in figure 9. The distribution of the molecular component of the scalar dissipation (
$\textit{Ra}|\boldsymbol{\nabla }C|^2$
) is shown in figures 9(a-i) and 9(a-ii), for a case without dispersion (
$\textit{Ra}=10^4,\ \varDelta \to \infty$
) and with dispersion (
$\textit{Ra}=10^4,\ \varDelta =10^{-1},\ r=10$
), respectively. The concentration fields refer to time
$t=10^4$
, at which the fingers have grown to less than half domain height (only a small region in the core of the domain is shown in figure 9
a,b). In the absence of dispersion (
$\varDelta \to \infty$
; figure 9
a-i), the fingers are constrained to grow vertically within the narrow space between two neighbouring fingers. This produces large values of dissipation at their side interface and also at their tips. Other regions of high dissipation are identified as the portion of domain near the centreline, where the finger pattern is more complex. In the dispersive case (
$\varDelta =10^{-1}$
; figure 9
a-ii), high values of dissipation are still localised in the same regions (side interface, tips and near the domain centreline), but corresponding values are much lower. This matter is further discussed from a global perspective in § 5.2.
A long-term consequence of the lower growth of the mixing region (
$\gamma =0.49$
) compared with the case without dispersion (
$\gamma =0.59$
) is that it will take longer for the fingers to reach the walls, and then to achieve the maximum value of mean scalar dissipation. When this happens, however, the subsequent dynamics is similar to that described in the unstably stratified and final diffusive regimes, leading to a progressive homogenisation of the concentration distribution. In the following section, the global mixing dynamics is presented.
5.2. Mixing
A behaviour similar to that discussed in § 4 is observed also in this case, with some variations: after an initial diffusion-dominated phase still well described by (4.3),
$\chi _m$
reported in figure 10(a) decreases when
$\varDelta$
is reduced, as a result of the thickening of the finger interface, leading to smoother gradients of concentration. The same trend was observed in previous studies in different flow configurations (Dhar et al. Reference Dhar, Meunier, Nadal and Meheust2022).
Evolution of the mean scalar dissipation for different values of
$\varDelta$
. The red line refers to the case in the absence of dispersion. The (a) molecular (
$\chi _m$
), (b) dispersive (
$\chi _d$
) and (c) total (
$\chi =\chi _m+\chi _d$
) dissipation. The initial diffusive solution (4.3) (dotted line) is also indicated.
Understanding the distribution of the dispersive scalar dissipation within the domain,
$\textit{Ra}[ (\boldsymbol{\nabla }C)\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C)-|\boldsymbol{\nabla }C|^2]$
reported in figure 9(b), is less trivial. In general, large values are found where both concentration gradient and velocity are large. In the absence of dispersion (figure 9
b-i), large (in magnitude, either positive or negative) values of dispersive dissipation correspond to finger tips and core. In contrast, when dispersion is considered (figure 9
b-ii), the whole mixing region is characterised by non-negligible values of dispersive dissipation. A more quantitative evaluation is provided by the probability density functions of the components of the molecular, dispersive and total dissipation relative to the specific fields considered, and reported in figure 9(c). Globally the results appear very different: without dispersion (figure 9
c-i) the distribution of dispersive dissipation is nearly symmetric and gives zero as mean global value, while in the dispersive case (figure 9
c-ii) the probability density function is skewed towards positive values. As a result, we have that the smaller the value of
$\varDelta$
, the larger the values of dispersive dissipation (see figure 10
b). Since the flow is initialised as still (
$\boldsymbol{u}=0$
), at early times the dispersion tensor (2.12) corresponds to
$\unicode{x1D63F}=\unicode{x1D644}$
, leading to a dispersive dissipation that is zero. As soon as convective instabilities form,
$\unicode{x1D63F}$
grows in magnitude, as well as
$\chi _d$
. The convective onset occurs earlier when dispersion is considered (see figure 10
b,c). In agreement with previous findings obtained in a different configuration (Dhar et al. Reference Dhar, Meunier, Nadal and Meheust2022), we find that for this large value of
$\textit{Ra}$
the effect of dispersion on the onset time is not very significant, while it may be much more pronounced at lower
$\textit{Ra}$
. In addition, the onset time is also non-monotonic with
$\varDelta$
. Eventually, for long times, the flow achieves a stably stratified configuration, and the strength of convective instabilities diminishes. As a result,
$\unicode{x1D63F}$
reduces again leading to a decrease of the total dissipation.
The total mean dissipation
$\chi$
is obtained combining the molecular (
$\chi _m$
, figure 10
a) and dispersive (
$\chi _d$
, figure 10
b) contributions, and it is reported in figure 10(c). Despite an earlier onset observed when dispersion is considered, which makes
$\chi$
depart sooner from the initial diffusive dissipation ((4.3), dotted line), simulations with large
$\varDelta$
initially exhibit larger values of total dissipation. For instance, it takes about
$t\approx 8\times 10^3$
for the case
$\varDelta =0.1$
to achieve values of
$\chi$
comparable to those at larger
$\varDelta$
, and even later the maximum total dissipation achieved is smaller compared with cases without dispersion. However, the situation changes dramatically after the fingers reach the walls, corresponding to the maximum of
$\chi$
in figure 10(c), between
$t\approx 2\times 10^4$
and
$t\approx 5\times 10^4$
, depending on
$\varDelta$
. For
$\varDelta =0.1$
and
$\varDelta =0.05$
, despite
$\chi$
being initially smaller compared with larger
$\varDelta$
, the total dissipation continues to contribute to mixing for a much longer time. Indeed, after the fingers reach the walls, the velocity reduces partially, but not suddenly. This provides a significant advantage compared with the cases without dispersion, especially when it comes to the cumulative amount of solute mixed, quantified by
$M$
.
Evolution of the degree of mixing (
$M$
) for different
$\varDelta$
and
$\textit{Ra}=10^4$
(
$r=10$
; see table 1 for further details). (a) Results are shown in terms of total mixing
$M$
, where a close-up view of the early phase is also reported in the inset. Here, the circle symbol marks the first instant considered in the simulations. The initial diffusive solution (4.3) (dotted line) is also indicated. The case without dispersion (
$\varDelta \to \infty$
, red line) is shown as a reference. (b) The relative importance of molecular mixing to total mixing,
$M_m/M$
, evaluated at each instant. Unsurprisingly, molecular mixing becomes progressively less important as
$\varDelta$
increases.
The degree of mixing
$M$
, reported in figure 11(a), confirms the previous observations. Initially all simulations considered follow the initial diffusive solution (4.3). Later, mixing in non-dispersive systems is more efficient than in cases without dispersion: the smaller the value of
$\varDelta$
, the sooner
$M$
will differ from the analytical diffusive case (inset of figure 11
a). The behaviour of
$M$
observed for
$\varDelta \geqslant 1$
is nearly independent of
$\varDelta$
for long times (
$t\geqslant 1\times 10^5$
; main panel of figure 11
a), attaining
$M\approx 0.85$
. When
$\varDelta$
is further diminished, namely
$\varDelta \leqslant 0.1$
, the evolution is impacted dramatically and the degree of mixing achieves values of 0.91 and 0.94 for
$\varDelta =0.05$
and
$\varDelta =0.01$
, respectively. This non-monotonic behaviour suggests that a non-trivial interplay of flow evolution and onset time controls the combined dynamics of molecular and dissipative mixing. On the one hand, the fraction of mixing due to molecular dissipation,
$M_m/M$
(figure 11
b), increases with
$\varDelta$
. This is expected, since the effect of reducing
$\varDelta$
is to increase the importance of dispersion. On the other hand, the flow configuration is such that the dispersive dissipation, produced by the interplay of velocity and concentration fields, is considerably larger for
$\varDelta =0.1$
than for
$\varDelta =0.05$
(see figure 10
b). This complex combination of different processes determines the evolution of the degree of mixing. Ultimately, all cases will achieve the well-mixed condition (
$M=1$
).
In conclusion, the presented results suggest that in this configuration (
$r=10,\ \textit{Ra}=10^4$
), within the range of
$\varDelta$
and times explored,
$\varDelta =0.1$
provides the most favourable mixing conditions and
$M$
achieves the largest values. The observation that there is an optimum value of
$\varDelta$
that maximises the mixing is opposite to previous findings in Rayleigh–Bénard configuration (Wen et al. Reference Wen, Chang and Hesse2018), where the flux was observed to be minimum for
$\varDelta =0.05$
. This is not surprising, given the very different nature (transient versus steady state, no external driving versus constant external driving) of the two systems considered.
6. Flow evolution with dispersion: influence of
$\boldsymbol{r}$
We analyse here the effect of the anisotropy ratio
$r$
, defined in (2.13), on the flow evolution and mixing dynamics. We fix
$\varDelta =0.1$
, corresponding to a strongly dispersive case studied in § 5, and we vary
$r$
in the interval
$1\leqslant r \leqslant 20$
. As previously mentioned,
$r$
varies with the flow parameters, and the most common practice to consider
$r\approx 10$
may not be motivated, unless the problem considered is solute transport in the advection-dominated regime (Bijeljic & Blunt Reference Bijeljic and Blunt2007).
Concentration fields at
$t=2\times 10^4$
for (a–f) different values of the dispersion parameter
$r$
. The field (a) corresponds to the case without dispersion (
$\varDelta \to \infty$
). See supplementary movies S2 and S3 for the time-dependent evolution of the simulations with
$r=10$
(e) and
$r=1$
(b).
6.1. Flow dynamics
The flow evolution follows the same regimes discussed in § 4 with some differences that are discussed in the following (see supplementary movies S2 and S3 for the time-dependent evolution of the simulations with
$r=10$
and
$r=1$
reported in figures 7(e) and 7(b), respectively). The concentration fields taken at time
$t=2\times 10^4$
are reported in figure 12(b–f) for
$\varDelta =0.1$
and
$1\leqslant r\leqslant 20$
, and compared with the case without dispersion (figure 12
a). We observe that, at this value of
$\varDelta$
, the anisotropy ratio
$r$
does have a role in shaping the finger pattern: the morphology of the flow changes dramatically when
$r$
is increased, corresponding to a coarsening of the fingers that become also more intricate near the centreline as
$r$
rises, in agreement with previous observations by Ghesmat & Azaiez (Reference Ghesmat and Azaiez2008).
Evolution of the horizontally averaged concentration profiles,
$\overline {C}$
, relative to simulation
$\textit{Ra}=10^4$
,
$\varDelta =0.1$
and
$r=1$
. Profiles reported correspond to instants taken preceding the finger impact on the walls. Specifically, they are in the range
$400\leqslant t\leqslant 16000$
(a), in the initial diffusive regime (
$t\leqslant 1.4\times 10^3$
) (b) and in the finger merging and growth regime (
$t\geqslant 7\times 10^3$
) (c). In (b), the dashed line indicates the initial diffusive solution (4.1). In (c), the wall-normal coordinate is rescaled with
$t-t_0$
, where
$t_0=4\times 10^3$
. The dashed line represents (4.4).
At early times (
$t\leqslant 1400$
), the flow dynamics is unaffected by dispersion since in this phase
$\boldsymbol{u}\approx 0$
. Therefore, the concentration field evolves according to the diffusive self-similar solution (4.1). This is confirmed by the horizontally averaged concentration profiles relative to the simulation with
$r=1$
and shown in figure 13(b). At later times, the fingers form and interact. The concentration fields in figure 12 suggest that the concentration gradients become weaker as
$r$
is increased.
We observe that, similarly to the case
$r=10$
(figure 8
c), the mixing layer grows approximately linearly. Indeed, for
$t\geqslant 7\times 10^3$
, the horizontally averaged concentration profiles of figure 13(c) are well fitted by the linear function (4.4) with
$t = 4\times 10^3$
and
$\gamma =0.46$
, where the growth rate of the mixing region,
$\gamma$
, matches the value obtained for
$r=10$
(discussed in § 5.1). This observation suggests that it will take a longer time for the fingers to reach the walls (corresponding to the instant at which the maximum of dissipation is achieved) compared with the case with
$\varDelta \to \infty$
. Following the finger’s impact on the horizontal boundaries, the flow evolves towards a progressive homogenisation of the concentration distribution.
Evolution of the mean scalar dissipation for different values of
$r$
. The red line refers to the case in the absence of dispersion. The (a) molecular (
$\chi _m$
), (b) dispersive (
$\chi _d$
) and (c) total (
$\chi =\chi _m+\chi _d$
) dissipation. The initial diffusive solution (4.3) (dotted line) is also indicated.
6.2. Mixing
At early times, the components of the scalar dissipation reported in figure 14 follow the purely diffusive solutions (4.3). We observe that the smaller the value of
$r$
, the sooner the onset of these instabilities occurs, making the dispersive curves in figure 14(b,c) deviate from the diffusive solutions. The onset is clearly driven by the dispersive component of the dispersion tensor: the molecular dissipation (figure 14
a) follows the diffusive solution for very long times, while the dispersive dissipation (figure 14
b) is now positive, with smaller values of
$r$
leading to an earlier growth compared with larger values of
$r$
, in agreement with previous findings (Dhar et al. Reference Dhar, Meunier, Nadal and Meheust2022). After formation, fingers grow vertically eventually interacting with neighbouring flow structures. The concentration gradients across the finger interface reduce as
$r$
is increased (figure 12), and so does the molecular component of the dissipation,
$\chi _m$
, as appears from figure 14(a).
As discussed in § 6.1, the fingers grow at an approximately constant rate (
$\gamma =0.46$
) that is lower compared with the case without dispersion (
$\gamma =0.59$
). We can infer this also from the dissipation curves in figure 14: the maximum of each curve is achieved at nearly the same time (
$t\approx 3\times 10^4$
). However, the maximum value of dispersive dissipation
$\chi _d$
is remarkably affected by
$r$
, as one would expect looking at the form of the dispersion tensor (2.12). For
$\varDelta =0.1$
considered here, when
$t\geqslant 4\times 10^3$
the dispersive component of the mean scalar dissipation is always dominant compared with the molecular counterpart. As a result of the interplay between molecular and dispersive dissipation, the behaviour of the total dissipation
$\chi$
, reported in figure 14(c), suggests that initially the mixing process is less efficient compared with the case without dispersion (
$\varDelta \to \infty$
, red line). In turn, for long times and after the fingers have reached the walls,
$\chi$
remains higher in the cases with dispersion, with important consequences on the degree of mixing.
Evolution of the degree of mixing (
$M$
) for different
$r$
and
$\textit{Ra}=10^4$
(
$\varDelta =0.1$
; see table 1 for further details). (a) Results are shown in terms of total mixing
$M$
, where a close-up view of the early phase is also reported in the inset. Here, the circle symbol marks the first instant considered in the simulations. The initial diffusive solution (4.3) (dotted line) is also indicated. The case without dispersion (
$\varDelta \to \infty$
, red line) is shown as a reference. (b) The relative importance of molecular mixing to total mixing,
$M_m/M$
, evaluated at each instant. Also in this case, molecular mixing becomes progressively less important as
$r$
increases.
The evolution of the degree of mixing (
$M$
) for different values of
$r$
is shown in figure 15. At very early times (
$t\leqslant 1500$
) (see the inset of figure 15
a), the degree of mixing is larger in the simulations with dispersion, due to the earlier onset of convection previously described. The fingers develop slower and within wider spaces compared with the case without dispersion: the gradients of concentration across the finger interface are smaller, corresponding to values of
$\chi _m$
that decrease with
$r$
. Later, the dispersive mixing increases in time and also with
$r$
. As a result of the interplay of these processes, the following dynamics appears: (i) the total degree of mixing is initially large for cases with dispersion (finite
$\varDelta$
) and low
$r$
; (ii) shortly after the fingers have formed,
$M$
becomes dominated by molecular diffusion and it is larger in the absence of dispersion (figure 15
a, main panel;
$2\times 10^3\leqslant t\leqslant 4.5\times 10^4$
); and (iii) after the fingers have reached the walls,
$\chi _d$
decreases at a lower rate compared with
$\chi _m$
, and
$M$
is controlled by the dispersive mixing. The relative importance of molecular mixing to total mixing,
$M_m/M$
, evaluated at each instant, is shown in figure 15(b). As expected, in the long term the molecular mixing becomes progressively less important as
$r$
is increased, namely from 46 % to 17 % for
$r=1$
and
$r=20$
, respectively.
7. Conclusions and outlook
7.1. Summary and conclusions
We analysed the process of convective mixing in the presence of dispersion in two-dimensional, homogeneous and isotropic porous media. We considered a Rayleigh–Taylor instability in which the presence of a solute produces density differences that drive the flow. The domain consists of two fluid layers, initially divided by a flat interface, which will mix driven by buoyancy forces and eventually reach a uniform solute distribution. Solute is redistributed by advection and dispersion. The effect of dispersion is modelled using the anisotropic Fickian dispersion tensor formulation (2.12) proposed by Bear (Reference Bear1961). In this model, in addition to molecular diffusion
$D_m^*$
, the solute is also redistributed by the spreading produced by the flow due to the complex pathways followed by the fluid parcels within the pore spaces (mechanical dispersion). This additional solute redistribution is quantified by the longitudinal (
$D_l^*$
, in the direction of the flow) and the transverse (
$D_t^*$
, perpendicular to the direction of the flow) dispersion coefficients. The flow is controlled by three dimensionless parameters: the Rayleigh–Darcy number
$\textit{Ra}$
(defining the relative strength of convection and diffusion) and the dispersion parameters
$r=D_l^*/D_t^*$
and
$\varDelta =D_m^*/D_t^*$
. With the aid of numerical Darcy simulations, we investigated the mixing dynamics without and with dispersion.
Influence of the dispersion parameters on the degree of mixing
$M$
:
$\varDelta$
(a) and
$r$
(b). Results are reported in terms of time degree of mixing
$M$
relative to the case without dispersion,
$M(\varDelta \to \infty )$
. The first instant considered in the simulations (circle) is also indicated.
We found that in the absence of dispersion (
$\varDelta \to \infty$
) the system’s dynamics is self-similar and independent of
$\textit{Ra}$
until the fingers approach the domain’s horizontal boundaries, and the flow evolution follows several regimes, which we describe. Then we analysed the effect of dispersion in time (
$t$
) for a fixed value of Rayleigh number (
$\textit{Ra}=10^4$
). We quantify the mixing state of the system using the degree of mixing
$M$
(Jha et al. Reference Jha, Cueto-Felgueroso and Juanes2011), which varies between
$M(t=0)=0$
(segregated layers) and
$M(t\to \infty )=1$
(uniform solute concentration field). A detailed analysis of the scalar dissipation reveals a complex interplay between flow structures and mixing due to the dispersive and molecular contributions. Both
$r$
and
$\varDelta$
affect the finger formation and development: in particular, the lower the value of
$\varDelta$
(or the larger the value of
$r$
), the wider, the more convoluted and diffused the fingers (see figures 7 and 12). For an anisotropic dispersion tensor with
$r=10$
(representative of solute transport in the advection-dominated regime; Bijeljic & Blunt Reference Bijeljic and Blunt2007), the role played by the relative importance of molecular and transverse dispersion,
$\varDelta$
, is crucial. This is presented in figure 16(a) in terms of degree of mixing,
$M(\varDelta )$
, normalised by the case without dispersion
$M(\varDelta \to \infty )$
. Three main phases appear: (i) initially (
$t\leqslant 2\times 10^3$
) the mixing is more efficient in cases with dispersion; (ii) after the fingers have grown sufficiently and merged, the mixing performance in the dispersive cases is lower than in the absence of dispersion; and (iii) for longer times, in contrast, high-dispersion flows (
$\varDelta \leqslant 0.1$
) show higher degree of mixing compared with systems without dispersion. Ultimately, all cases will lead to the same uniform configuration (
$M=1$
). Remarkably, we found that within the time frame considered, there is an optimum value
$\varDelta =0.1$
that maximises the mixing. We finally looked at the impact of
$r$
on the mixing when
$\varDelta =0.1$
(figure 16
b). We found that, for such value of
$\varDelta$
, the anisotropy ratio
$r$
produces only second-order effects, with some noticeable changes only in the intermediate phase, where it appears that the mixing is more efficient for small values of anisotropy.
7.2. Implications for geophysical flows: saline seepage in groundwater systems
The theoretical framework proposed allowed us to analyse numerical Darcy simulations with dispersion by splitting the mean scalar dissipation into molecular and dispersive components. In combination with pore-scale simulations and bead packs experiments, this approach can be used to validate and improve current dispersion models to obtain more reliable estimates of solute transport and spreading in the subsurface. As an example of a possible application, we discuss here the case of saline seepage from salt water basins presented by Narayan & Armstrong (Reference Narayan and Armstrong1995), and summarised in the following.
(a) Conceptualised hydrogeology of the River Murray basin area, adapted from Narayan & Armstrong (Reference Narayan and Armstrong1995). (b) Modelling of the saline seepage through the bottom of Lake Ranfurly West. (i) The high-permeability sands aquifer is confined by two low-permeability layers. (ii) Lake Ranfurly West supplies high-salt-concentration water (
$C^*=C^*_{\textit{max}}$
) from the top, while low-salinity (
$C^*=C^*_{\textit{min}}$
) groundwater is present in the aquifer.
One of the most important drainage systems in Australia is represented by the Murray-Darling river, a key source of water in the region. Near-surface groundwater in this basin is characterised by high salinity. Agricultural activities have led to a rising of the water tables, with the consequence of an increased discharge of high-salt-concentration groundwater into the river basin. This process may eventually increase the river’s water salinity to unacceptable levels in periods of low flow rate in the river. To prevent this issue, high-salinity groundwater is intercepted and stored in basins at the surface level, where it may evaporate, further increasing the salt concentration. In some cases, these surface basins are designed to allow a slow and controlled leak to the underlying Murray River aquifer. This is the case for Lake Ranfurly West (Ghassemi, Thomas & Jakeman Reference Ghassemi, Thomas and Jakeman1988), which releases high-salinity water through the Channel sands aquifer into the River Murray. The hydrogeology of the system is sketched in figure 17(a). Designing and controlling such basins are key to manage the water resources efficiently and to keep the salinity of the rivers at an acceptable level. Here we apply our findings to determine the role of dispersion in the salt spreading process from Lake Ranfurly West and the River Murray basin.
The flow can be modelled as a high-permeability region (Channel sands aquifer) confined by two low-permeability clay layers (see figure 17
b-i): the clay layer at the bottom of the lake confines the system from above, and the Upper Parilla clays confine the system from below (Narayan & Armstrong Reference Narayan and Armstrong1995). The fluid saturating the high-permeability region can be modelled as a heavy and high-salt-concentration fluid (
$C^*=C^*_{\textit{max}}=0.120\,{\rm kg\,kg}^{-1}$
, seeping from Lake Ranfurly West) sitting on top of a lighter, low-salinity fluid (
$C^*=C^*_{\textit{min}}=0.045\,{\rm kg\,kg}^{-1}$
, groundwater present in the Channel sands aquifer). A simplified flow configuration associated with this system is illustrated is figure 17(b). We neglect here any flow circulations induced by pumping, and we assume that the flow is only driven by density differences. We take the system parameters from Narayan & Armstrong (Reference Narayan and Armstrong1995) and determine the relevant dimensionless flow quantities as follows. The depth of the porous layer is
$L_z^*=4$
m, with porosity
$\phi =0.3$
and permeability
$k=2.95\times 10^{-11}$
m
$^2$
(we consider the medium isotropic). The salt concentration difference
$C^*_{\textit{max}}-C^*_{\textit{min}}$
produces a maximum density difference
$\Delta \rho ^*=52.5$
kg m
$^{-3}$
. Considering a dynamic viscosity
$\mu =10^{-3}$
Pa s, the characteristic velocity is
$\mathcal{U}^*=1.52\times 10^{-5}$
m s–1
$=1.31$
m d–1, and we obtain a Rayleigh number
$\textit{Ra}=1.35\times 10^5$
. The longitudinal dispersivity is
$\alpha ^*_l = 80$
m and the anisotropy ratio
$r=10$
. With this set of parameters, we have
$\varDelta = rD_m/(\mathcal{U}^*\alpha ^*_l)\approx 10^{-5}$
. From the results presented in figure 16(a) we conclude that dispersion is the dominant mixing mechanism, and has to be accounted for in the design and simulation of these flows. We also remark that, despite the high uncertainty on
$\alpha ^*_l$
, also much smaller values would lead to the same conclusion, e.g. for
$\alpha ^*_l = 0.8$
m we have that
$\varDelta \approx 10^{-3}$
, still in the range of parameters in which dispersion represents the dominant mixing mechanism. This would change only if substantially lower values of
$\textit{Ra}$
are considered.
7.3. Limitations and future developments
A possible limitation of the present work consists of the model adopted. Despite being widely employed (De Paoli Reference De Paoli2023), the dispersion model considered is derived for idealised advective flows, i.e. where the source of driving is a pressure gradient that generates a constant and uniform velocity field. However, the configuration studied here is transient, and therefore the flow velocity scales also evolve across the different regimes experienced by the flow. An interesting possible development consists of taking the variation of the velocity scales into account in the dispersion model, to better describe the dispersion coefficients
$D_l^*$
and
$D_t^*$
(Koch & Brady Reference Koch and Brady1985; Delgado Reference Delgado2007). In addition, recently a theory that quantifies the interplay between intrapore and interpore flow variabilities and their impact on hydrodynamic dispersion has been proposed (Liu et al. Reference Liu, Xiao, Aquino, Dentz and Wang2024). Extending the model of Liu et al. (Reference Liu, Xiao, Aquino, Dentz and Wang2024) to explore buoyancy-driven mixing would improve the reliability of the results obtained via numerical simulations in the present context (Woods Reference Woods2025).
The dimensionality of the system will likely have an effect on the mixing process, as previously observed by Boffetta et al. (Reference Boffetta, Borgnino and Musacchio2020) and Boffetta & Musacchio (Reference Boffetta and Musacchio2022), who studied the present flow configuration without dispersion. In particular, they have shown that the growth of the mixing layer is strongly affected by the dimensionality of the domain, with the fingers being more coherent in two than in three dimensions. This leads to a faster growth of the mixing region and to a larger variance of the concentration distribution in the two-dimensional case. The effects of dimensionality have been recently explored in the one-sided configuration (De Paoli et al. Reference De Paoli2025b
), and the same flow dynamics has been observed in two and in three dimensions, with a larger finger size and a lower mixing rate in two dimensions. In contrast, still in the one-sided configuration but in the presence of dispersion, Dhar et al. (Reference Dhar, Meunier, Nadal and Meheust2022) observed that the finger size is independent of the dimensionality at
$\textit{Ra}=10^3$
, while it is larger in three than in two dimensions at
$\textit{Ra}=3\times 10^3$
, and in all cases the fingers grow at the same rate. They also observed that the effect of dispersion on the onset time and on the mixing are similar in two and in three dimensions. Considering these findings, we expect that in a three-dimensional, porous and dispersive Rayleigh–Taylor flow, the mixing dynamics will exhibit: (i) a slower growth rate of the mixing region, (ii) wider fingers and (iii) similar effects of dispersion on the mixing dynamics, compared with the two-dimensional case presented here. To confirm these assumptions, a full campaign of three-dimensional numerical simulations is required.
In this work, we considered fully miscible fluids. In geological processes, however, the presence of partially miscible fluids (e.g. CO
$_2$
and brine) may significantly affect the flow dynamics and solute mixing processes. With the aid of a phase-field method (Darcy–Cahn–Hilliard), Li et al. (Reference Li, Cai, Chen and Meiburg2022, Reference Li, Lin, Cai, Chen and Meiburg2023) have recently shown that, in semi-infinite domains, partial miscibility of the fluids leads to larger dissolution fluxes compared with the case of fully miscible fluids, corresponding to a more efficient mixing process. We believe that in the present configuration the same effect could possibly occur at early times, when the driving, provided by the volume of fresh fluid available outside the mixing region, is nearly constant.
Finally, our findings have been obtained in idealised conditions, since we considered homogeneous and isotropic porous domains. However, geological formations are heterogeneous at multiple scales, requiring Darcy simulations to include non-uniformities in the permeability and porosity fields (Woods Reference Woods2015). An interaction between the scales of the heterogeneity and of the instability that controls the flow evolution has been observed (Benhammadi, Meunier & Hidalgo Reference Benhammadi, Meunier and Hidalgo2025). Heterogeneities can also influence the form of the dispersion tensor, and possibly lead to non-Fickian effects (Cala & Greenkorn Reference Cala and Greenkorn1986). A systematic study of Rayleigh–Taylor flows in heterogeneous porous domains at the Darcy scale, combined with pore-scale simulations and experiments, is advisable to shed new light on the role of heterogeneities, and eventually to improve and validate existing dispersion models.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10684.
Acknowledgements
B. Wen is gratefully acknowledged for providing the data of Wen et al. (Reference Wen, Chang and Hesse2018), used to validate the code with dispersion. J. Wood, Y. Yang and C. Hu are also acknowledged for useful discussions. Three anonymous referees are also acknowledged for their insightful comments and suggestions.
Funding
This project has received funding from the European Union’s Horizon Europe research and innovation programme under Marie Sklodowska-Curie grant agreement MEDIA no. 101062123. We acknowledge the EuroHPC Joint Undertaking for awarding the projects EHPC-BEN-2024B08-060 and EHPC-EXT-2024E02-122 to access the EuroHPC supercomputer MareNostrum5 hosted at the Barcelona Supercomputing Center (Spain).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Additional details of the dispersion tensor
In a three-dimensional reference frame, the dispersion tensor (2.12) with all the terms written explicitly reads (Zheng, Wang & Reference Zheng and Wang1999)
\begin{equation} \unicode{x1D63F}= \begin{bmatrix} 1+\frac {1}{\varDelta }\left [\lvert \boldsymbol{u}\rvert +(r-1)\frac {u^{2}}{\lvert \boldsymbol{u}\rvert }\right ] & (r-1)\frac {uv}{\varDelta \lvert \boldsymbol{u}\rvert } & (r-1)\frac {uw}{\varDelta \lvert \boldsymbol{u}\rvert } \\ (r-1)\frac {uv}{\varDelta \lvert \boldsymbol{u}\rvert } & 1+\frac {1}{\varDelta }\left [\lvert \boldsymbol{u}\rvert +(r-1)\frac {v^{2}}{\lvert \boldsymbol{u}\rvert }\right ] & (r-1)\frac {vw}{\varDelta \lvert \boldsymbol{u}\rvert } \\ (r-1)\frac {uw}{\varDelta \lvert \boldsymbol{u}\rvert } & (r-1)\frac {vw}{\varDelta \lvert \boldsymbol{u}\rvert } & 1+\frac {1}{\varDelta }\left [\lvert \boldsymbol{u}\rvert +(r-1)\frac {w^{2}}{\lvert \boldsymbol{u}\rvert }\right ] \\ \end{bmatrix}\!. \end{equation}
Note that on expanding the divergence term of the right-hand side of (2.8) we obtain
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C) &= \frac {\partial }{\partial x_{i}}\left ( D_{ij}\frac {\partial C}{\partial x_{j}}\right )= \frac {\partial D_{ij}}{\partial x_{i}}\frac {\partial C}{\partial x_{j}}+ D_{ij}\frac {\partial }{\partial x_{i}}\left (\frac {\partial C}{\partial x_{j}}\right ) \end{align}
\begin{align} \begin{split} &= D_{xx}\frac {\partial ^{2} C}{\partial x^{2}}+D_{yy}\frac {\partial ^{2} C}{\partial y^{2}}+D_{zz}\frac {\partial ^{2} C}{\partial z^{2}}+ \frac {\partial C}{\partial x} \left (\frac {\partial D_{xx}}{\partial x}+\frac {\partial D_{xy}}{\partial y}+\frac {\partial D_{xz}}{\partial z}\right )\\ &\quad + \frac {\partial C}{\partial y} \left (\frac {\partial D_{xy}}{\partial x}+\frac {\partial D_{yy}}{\partial y}+\frac {\partial D_{yz}}{\partial z}\right )+ \frac {\partial C}{\partial z} \left (\frac {\partial D_{xz}}{\partial x}+\frac {\partial D_{yz}}{\partial y}+\frac {\partial D_{zz}}{\partial z}\right )\\ &\quad +2\left ( D_{xy}\frac {\partial ^{2} C}{\partial x\partial y} +D_{yz}\frac {\partial ^{2} C}{\partial y\partial z} +D_{xz}\frac {\partial ^{2} C}{\partial x\partial z}\right)\!. \end{split} \end{align}
The right-hand side of (A3) reduces to the Laplacian form
$D_m^*{\nabla} ^{2}C$
when
$D_{ii}=D_m^*$
(with
$D_m^*$
uniform in space) and
$D_{ij, i\ne j}=0$
. If the dispersion tensor (2.12) is isotropic (
$r=1$
; see § 6), then (A3) takes the form
\begin{align} \begin{split} \boldsymbol{\nabla }\boldsymbol{\cdot }(\unicode{x1D63F}\boldsymbol{\nabla }C) &= \frac {\partial }{\partial x_i}\left (D_{ii}\frac {\partial C}{\partial x_i}\right )\\ &=\frac {\partial C}{\partial x} \frac {\partial D_{xx}}{\partial x}+ \frac {\partial C}{\partial y} \frac {\partial D_{yy}}{\partial y}+ \frac {\partial C}{\partial z} \frac {\partial D_{zz}}{\partial z}+D_{xx}\frac {\partial ^{2} C}{\partial x^{2}}+D_{yy}\frac {\partial ^{2} C}{\partial y^{2}}+D_{zz}\frac {\partial ^{2} C}{\partial z^{2}}, \end{split} \end{align}
where the convention for summation with respect to repeated indices is used. In the case where the coefficients
$D_{ii}$
are constant, the Laplacian form is again recovered.
Appendix B. Numerical details
B.1. Numerical treatment of the dispersive term
A possible discretisation of (2.8) consists of splitting the terms as follows:
\begin{equation} \underbrace {\left (\frac {\partial }{\partial t} - {\nabla} ^2\right )C}_{\text{Implicit terms}} = \underbrace {\boldsymbol{\nabla }\boldsymbol{\cdot }\left[\left (\unicode{x1D63F}-\unicode{x1D644}\kern1.5pt\right) \boldsymbol{\nabla }C\right ] - \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }C}_{\text{Explicit terms}}. \end{equation}
In this method, the term on the left-hand side of the equation is treated implicitly and the term on the right-hand side is treated explicitly (the explicit treatment of the dispersive term was also employed by Wen et al. (Reference Wen, Chang and Hesse2018)). However, this approach is very restrictive on the size of the time step. To ensure numerical stability, the time step must scale with the square of the minimum grid spacing. For finely resolved boundary layers, this makes the computation unrealistically expensive. In order to remedy this issue, the following strategy is employed: (i) first, a predictor–corrector type of approach is used with (B1) serving as the predictor step, in order to obtain the intermediate solution for the concentration field
$C^{n+1/2}$
, with the superscripts
$n$
indicating current time step; (ii) then (2.9) and (2.10) are solved to obtain an intermediate solution for the velocity field
$\boldsymbol{u}^{n+1/2}$
; (iii) the dispersion coefficients for
$\unicode{x1D63F}^{\,n+1/2}-\unicode{x1D644}$
are computed using the intermediate velocity field; and (iv) the corrector step is performed by treating only the advective terms
$\boldsymbol{u}^n\boldsymbol{\cdot }\boldsymbol{\nabla }C^n$
explicitly, while the rest of the terms can be treated implicitly as
\begin{equation} \underbrace {\left [\frac {\partial }{\partial t} - \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\unicode{x1D63F}\boldsymbol{\nabla }\right )\right ]C}_{\text{Implicit terms}} = \underbrace { - \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }C}_{\text{Explicit terms}}. \end{equation}
Traditionally, the discretised form of (B1) could be solved in AFiD (De Paoli et al. Reference De Paoli2025a
) by factorising the implicit terms effectively splitting the corresponding discrete operator into three discrete linear operators along each of the Cartesian dimensions. Such an algorithm cannot be used to solve the discretised form of (B2), and therefore we use iterative methods. Since it is difficult to realise highly scalable parallel implementations of methods with fast convergence (e.g. Gauss–Seidel), the Jacobi iteration method (Ferziger, Perić & Street Reference Ferziger, Perić and Street2019) is used. For sufficiently large
$\textit{Ra}$
, the discretised operator matrix resulting from the implicit terms is strictly diagonally dominant, thereby guaranteeing convergence. Implicitly treating the dispersion terms in this way ensures the numerical stability of the solution for larger time steps. In this case, the size of the time step is only restricted by the Courant–Friedrichs–Lewy number determined by the advective terms.
In the following, the discretisation of the dispersive terms is explained by using a two-dimensional example in
$x,y$
space, but it can easily be extended to three dimensions. Consider a staggered grid as shown in figure 18 with grid points containing the information of the concentration field
$C$
at the cell centres. We use linear interpolations of the velocity components to compute
$D_{xx}$
at face centres where
$u$
is stored,
$D_{xy}$
at cell vertices and
$D_{yy}$
at face centres where
$v$
is stored (additional details of the variable arrangement on the grid are provided by De Paoli et al. (Reference De Paoli2025a
)). For the predictor step,
$D_{xx} - 1$
is used instead of
$D_{xx}$
,
$D_{yy} - 1$
is used instead of
$D_{yy}$
and
$D_{zz} - 1$
is used instead of
$D_{zz}$
. For the corrector step, we compute and store the coefficients of the 26 neighbouring concentration grid points in a stencil that spans three grid cells along each dimension. These coefficients are used as matrix elements in the Jacobi iteration method.
Grid layout for demonstrating the computation of dispersion terms in two dimensions. Grid points containing information of variable
$C$
are indicated in blue,
$D_{xx}$
in green,
$D_{yy}$
in red and
$D_{xy}$
in yellow. Note that the subscripts
$i,j$
here no longer refer to the indices of the general dispersion tensor
$\unicode{x1D63F}$
, instead indicating the grid coordinates along
$x$
and
$y$
axes, respectively. Additional details of the variable arrangement on the grid are provided by De Paoli et al. (Reference De Paoli, Yerragolam, Lohse and Verzicco2025a
).
B.2. Initial condition and grid resolution
The initial condition is determined from the analytical solution (4.1) of the advection–dispersion equation (2.8), discussed in § 4. In correspondence of the time (
$t_0$
) employed to initialise the flow, the initial concentration distribution is
\begin{equation} C(z,t=t_0)=\frac {1}{2}\left [1+\text{erf}\left (\frac {z}{2\sqrt {t_0}}\right )\right ]\!. \end{equation}
The inverse function associating
$\tilde {z}$
to a given value of concentration
$\tilde {C}$
at time
$t_0$
is
\begin{equation} \tilde {z}(\tilde {C},t=t_0)=2\sqrt {t_0}\text{erf}^{\,-1}(2\tilde {C}-1) \end{equation}
and allows one to determine the minimum requirements in terms of combination of grid resolution/initial time considered. Assuming we aim at solving from 1 % to 99 % of the initial scalar difference with at least
$n_z$
points and the grid used is uniform, using (B4) we have that the grid spacing in the wall-normal direction is
\begin{equation} \Delta z = \frac {\textit{Ra}}{N_z-1} \leqslant \frac {\tilde {z}_{99}-\tilde {z}_1}{n_z} = \frac {2\sqrt { t_0}}{n_z}\left [\text{erf}^{\,-1}(0.98)-\text{erf}^{\,-1}(-0.98)\right ] \approx 6.580\frac {\sqrt {t_0}}{n_z}. \end{equation}
Therefore, to resolve the 1 % to 99 % scalar variation with at least
$n_z$
grid points we need to initialise the flow with
$t_0\geqslant 0.0231\times [n_z\textit{Ra}/(N_z-1)]^2$
. Using
$N_z=1024$
and
$\textit{Ra}=10^4$
, it gives
$t_0\approx 20$
(
$n_z=3$
) and
$t_0\approx 80$
(
$n_z=6$
). All simulations are initialised with the same initial perturbation at
$t_0=50$
. The simulations are initially run without dispersion up to
$t=200$
, and the dispersive terms are later activated. Note that
$t=200$
falls well within the diffusive part of the simulations (see e.g. figure 10
a), and therefore this formulation does not affect the subsequent formation and development of the flow structures. In addition, this strategy ensures that the same initial concentration and velocity fields are employed in all the simulations, making the different behaviour observed for
$r$
and
$\varDelta$
independent of the initial flow configuration, and also providing a smooth initial condition to the simulations with dispersion. We remark that minor discrepancies in the diffusive regime compared with the theoretical solution (see
$\chi _d$
at early times and for high dispersion; figures 10
b and 14
b) are due to the resolution, and doubling the resolution would resolve this minor issue.
Appendix C. Analytical solution for diffusion in confined domains
C.1. Initial diffusion at low Rayleigh number
We derive here the analytical solution for the mean scalar dissipation assuming pure diffusion in a confined domain with no-flux boundary conditions. We use (2.8), we assume that the fluid is still (
$\boldsymbol{u}=0$
) and we consider the system uniform in the horizontal direction (
$\partial _x C = 0$
). In addition, we assume that this regime starts at
$t=t_i$
. Introducing the rescaled dimensionless coordinate
$\bar {z} = z + \textit{Ra}/2$
, the exact solution is (Strauss Reference Strauss2007)
\begin{equation} C(\bar {z},t) = \sum _{n=0}^\infty A_n \exp {\left [-\left (\frac {n\pi }{\textit{Ra}}\right )^2t\right ]}\cos {\left (\frac {n\pi \bar {z}}{\textit{Ra}}\right )}\!\text{ .} \end{equation}
The coefficients of the cosine series (C1) are defined as
\begin{equation} A_n = \begin{cases} \displaystyle \frac {1}{\textit{Ra}}\int _0^{\textit{Ra}} C(\bar {z},t_i) \textrm { d}\bar {z},\quad \text{for } n=0, \\[10pt] \displaystyle \frac {2}{\textit{Ra}}\int _0^{\textit{Ra}} C(\bar {z},t_i) \cos {(n\pi \bar {z})} \textrm { d}\bar {z},\quad \text{for } n\geqslant 1, \end{cases} \end{equation}
with
$C(\bar {z},t_i)$
defined as the initial condition (B3), which in terms of
$\bar {z}$
reads
\begin{equation} C(\bar {z},t_i) = \frac {1}{2}\left [1+\text{erf}\left (\frac {\bar {z}-\textit{Ra}/2}{2\sqrt {t_i}}\right )\right ]. \end{equation}
Using the definition of molecular mean scalar dissipation (3.10) and the analytical diffusive solution (C1), the mean scalar dissipation is
\begin{align} \chi _m(t+t_i)&=\int _{0}^{\textit{Ra}}\left (\frac {\partial C}{\partial \bar {z}}\right )^2 \text{ d}\bar {z} \end{align}
\begin{align} &=\int _0^{\textit{Ra}}\left (\sum _{n=0}^\infty \frac {n\pi }{\textit{Ra}} A_n \exp {\left [-\left (\frac {n\pi }{\textit{Ra}}\right )^2t\right ]}\sin {\left (\frac {n\pi \bar {z}}{\textit{Ra}}\right )}\right )^2\text{ d}\bar {z}. \end{align}
The solution (C5) for
$\textit{Ra}=10^2$
is shown in figure 2. It has been computed using
$n=100$
, with a spatial discretisation of 128 points for
$\bar {z}$
and
$t_i=50$
, i.e. the same grid size and initial time used in the corresponding simulation. Note the time shift
$t_i$
in the solution (C5), due to the fact that (C1) is computed assuming
$t=t_i$
as the initial time. The accuracy of the analytical solution in predicting the decay of dissipation obtained in the simulation is excellent. In addition, we observe in (C5) that
$\partial C/\partial \bar {z}\sim \exp (-n^2)$
, indicating a fast decay with
$n$
. We verified that
$n=1$
is sufficient to capture well the decay, and additional terms improve the behaviour only for the very early times
$(t\leqslant 60)$
.
C.2. Final diffusion at high Rayleigh numbers
As discussed in § 4, the evolution during the final diffusive regime starts at time
$t_f$
with a linear concentration profile:
\begin{equation} C(\bar {z},t_f) = \frac {1}{2} + \beta \left (\frac {1}{2} - \frac {\bar {z}}{\textit{Ra}}\right )\!. \end{equation}
The initial condition is set with
$\beta =0.25$
and
$t_f=2.5\times 10^4$
for
$\textit{Ra}=10^3$
and
$\beta =0.70$
and
$t_f=10^5$
for
$\textit{Ra}=10^4$
, and it is used to determine the cosine series coefficients as in (C2). Finally, the evolution of the mean scalar dissipation is determined using (C5) as
\begin{equation} \chi _m(t+t_f) = \int _0^{\textit{Ra}}\left (\sum _{n=0}^\infty \frac {n\pi }{\textit{Ra}} A_n \exp {\left [-\left (\frac {n\pi }{\textit{Ra}}\right )^2t\right ]}\sin {\left (\frac {n\pi \bar {z}}{\textit{Ra}}\right )}\right )^2\text{ d}\bar {z}, \end{equation}
and the results are reported in figure 6, where a spatial discretisation of 512 points for
$\bar {z}$
is used.
Open access
