2. Mathematical description of flow and transport

2.1. Governing equations

In this work we consider a two-dimensional porous medium in a dimensionless Cartesian coordinate system ${\boldsymbol {x}} = (x, z)$. The dimensionless lengths of our domain are $L$ in the horizontal direction and $H$ in the vertical direction. All equations and quantities in this section are normalised and the dimensionless groups employed are reported in Appendix A. In particular, lengths are normalised by a characteristic length $\ell$ of the domain, chosen here to be the thickness of the domain. Thus, we have $H=1$. We start by considering a homogeneous porous medium, characterised by a constant permeability $\hat {k}$. Assuming incompressible flow, in the absence of temporally variable boundary conditions and sources and sink terms, the governing (dimensionless) equation for the flow field is given by

(2.1)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q}(\boldsymbol{x},t) = 0, \end{equation}

where $\boldsymbol {q}$ denotes the specific discharge, defined by Darcy's law as

(2.2)\begin{equation} \boldsymbol{q}(\boldsymbol{x},t) ={-} \frac{1}{\mu(c(\boldsymbol{x},t))} [\boldsymbol{\nabla} p(\boldsymbol{x},t) - \mathcal{F}(c(\boldsymbol{x},t))\boldsymbol{e}_z ]. \end{equation}

The velocity field is given by $\boldsymbol {q}/\phi$, where $\phi$ is the medium porosity. In (2.2), $p({\boldsymbol {x}},t)$ denotes the pressure field and $\mu (c)$ the viscosity of the ambient fluid mixture of concentration $c$. The unit vector along the vertical $z$ direction is denoted by $\boldsymbol {e}_z$. For our work, we adopt an exponential viscosity model (Jha et al. Reference Jha, Cueto-Felgueroso and Juanes2011a) to relate $\mu$ and $c$. An exponential viscosity model can capture the concentration dependence of viscosity of mixtures such as water and glycerol (Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996). The dimensionless form of the functional relationship between $\mu$ and $c$ used in this work is

(2.3)\begin{equation} \mu(\boldsymbol{x},t) = \mathrm{exp}[R(1-c(\boldsymbol{x},t))], \end{equation}

where $R$ is the log-viscosity ratio $R=\mathrm {ln}(\mu _0/\mu _1)$, in which $\mu _0$ is the viscosity of the ambient fluid and $\mu _1$ is the viscosity of the injected solute. For reference, in geological carbon sequestration applications, the log-viscosity ratio $R$ ranges between 1.9 and 2.7 for the viscosity values of CO$_2$ and water reported in Frailey & Leetaru (Reference Frailey and Leetaru2009).

The function $\mathcal {F}(c)$ in (2.2) originates from the density profile model adopted in Daniel & Riaz (Reference Daniel and Riaz2014):

(2.4)\begin{equation} \rho(c(\boldsymbol{x},t)) = \rho_0 + \Delta \rho \mathcal{F}(c(\boldsymbol{x},t)), \end{equation}

where $\rho _0$ is the ambient fluid density (i.e. the density at $c=0$) and $\Delta \rho = \rho _m -\rho _0$, in which $\rho _m$ is the maximum possible density of the mixture. In our work, we set $F(c)=c$, thus adopting a linear density profile (see Horton & Rogers Reference Horton and Rogers1945; Daniel & Riaz Reference Daniel and Riaz2014; Cowell, Kent & Trevelyan Reference Cowell, Kent and Trevelyan2020). A linear density profile is able to describe the density dependence on concentration of a potassium permanganate and water mixture; such a mixture has been employed in experimental studies such as that of Slim et al. (Reference Slim, Bandi, Miller and Mahadevan2013) as a model for $\mathrm {CO}_{2}$ in brine. For the linear model employed, the maximum density $\rho _m$ will coincide with the density of the injected solute $\rho _1$ (i.e. the density at $c=1$). Therefore (2.2) becomes

(2.5)\begin{equation} \boldsymbol{q}(\boldsymbol{x},t) ={-} \frac{1}{\mu(c(\boldsymbol{x},t))} [\boldsymbol{\nabla} p(\boldsymbol{x},t) - c(\boldsymbol{x},t)\boldsymbol{e}_z ]. \end{equation}

The governing dimensionless transport equation for an inert substance is given by

(2.6)\begin{equation} \frac{\partial c(\boldsymbol{x},t)}{\partial t} + \boldsymbol{q}(\boldsymbol{x},t) \boldsymbol{\cdot} \boldsymbol{\nabla} c(\boldsymbol{x},t) - \frac{1}{Ra} \nabla^2 c(\boldsymbol{x},t) =0, \end{equation}

in which the Rayleigh number $Ra$ is defined as

(2.7)\begin{equation} Ra = \frac{k_c \Delta \rho g \ell}{\phi D \mu_1}, \end{equation}

where $k_c$ is the medium characteristic permeability (for a homogeneous medium, $k_c=\hat {k}$), $\ell$ is the domain characteristic length (here, $\ell$ represents the vertical dimension of the porous formation), $\phi$ is the medium porosity and $D$ is the isotropic local-scale dispersion coefficient. Both $\phi$ and $D$ are assumed to be constant.

2.2. Boundary and initial conditions

The boundary conditions for the flow field are, as shown in figure 1, a no-flux condition on the top and bottom of the domain, a constant inlet flux $f_L$ on the left-hand boundary and an assigned constant pressure $p_R$ on the right-hand boundary. These are mathematically equivalent to

(2.8)$$\begin{gather} \left.\frac{\partial p}{\partial z}\right|_{z=0}=\left.\frac{\partial p}{\partial z}\right|_{z=H}=0, \end{gather}$$

(2.9)$$\begin{gather}\left.\frac{\partial p}{\partial x}\right|_{x=0}=f_L, \end{gather}$$

(2.10)$$\begin{gather}p|_{x=L}=p_R. \end{gather}$$

The dimensionless inlet flux $f_L$ can be understood as a measure of the relative strength of horizontal advection compared with diffusion (equivalent to a Péclet number), whereas $Ra$ controls the relative strength of vertical convection compared with diffusion. Therefore, we expect both parameters to control the strength of convective instabilities – viscous fingering and gravitational fingering – with the effect of $Ra$ being limited to the latter.

For the transport problem, we adopt a no-flux boundary condition at the top and bottom edges of the computational domain, analogously to the flow problem, which is consistent with the assumption of confining impermeable layers:

(2.11)\begin{equation} \left.\frac{\partial c}{\partial z}\right|_{z=0}=\left.\frac{\partial c}{\partial z}\right|_{z=H}=0. \end{equation}

The no flux-boundary condition applies to the left-hand boundary as well:

(2.12)\begin{equation} \left.\frac{\partial c}{\partial x}\right|_{x=0}=0, \end{equation}

while on the right-hand boundary of the domain we assume a natural outflow boundary condition, i.e. the solute mixture is allowed to exit the domain through convection only, following the direction of the average flow velocity.

The initial condition for the concentration field of the injected fluid is a zero concentration in the whole domain except where the initial source is located. Our initial source zone is represented by an ellipse of radii $a$ and $b$ along the longitudinal and transverse directions respectively, as displayed in figure 1. Within the source zone, represented here by the area $\mathcal {S}$ of the ellipse, the value of concentration is unitary at initial time $t=t_0$:

(2.13)\begin{equation} c(\boldsymbol{x},t_0) = \begin{cases} 1, & \mathrm{for}\ \boldsymbol{x} \in \mathcal{S},\\ 0, & \mathrm{for}\ \boldsymbol{x} \notin \mathcal{S}. \end{cases} \end{equation}

By changing the ratio $a/b$, we can explore the effects of source orientation (figure 1) on the spatiotemporal dynamics of $c$ and mixing metrics. Note that the area of the ellipse $\mathcal {S}$ is maintained constant throughout all the values of the ratio $a/b$ explored.

The theory of linear stability analysis of viscosity- and density-driven fingering can provide the instability onset time and the wavenumber versus growth rate relation for the early-time behaviour of the system. The linear analysis, including the non-modal analysis, is valid before the fingers start to interact via nonlinear mechanisms such as tip-splitting, shielding, coalescence, coarsening and channelling, which control plume spreading and mixing. Given our focus on spreading and mixing, we use direct numerical simulations, instead of the linear theory, to resolve such nonlinear interactions and their impact on spreading and mixing metrics. It is important to note that in our numerical model configuration (figure 1), the background flow direction is not parallel with the direction where the effects of gravity are manifested. To the best of the authors’ knowledge, the literature does not provide results from the linear stability analysis that allows for the prediction of instability wave number and growth rate for this specific configuration. Due to the orthogonality of background flow and gravity direction in our configuration, macroscopic vorticity can arise, as visible in figure 2. It is possible to observe in figure 2 that, in the case of a horizontal initial source characterised with multiple density fingers, the faster-moving left-side finger grows downward to reach the bottom boundary while the slower-moving right-side finger grows leftward to merge with the rest of plume at the lower boundary. As a result, the plume appears to roll over itself like a vortex. Merging of viscous fingers is a well-known mechanism (Jha et al. Reference Jha, Cueto-Felgueroso and Juanes2011b, Reference Jha, Cueto-Felgueroso and Juanes2013), but rolling of fingers, as observed in figure 2, due to the interplay between gravitational and viscous fingering mechanisms has not been mentioned in the literature.

Figure 2.

Temporal evolution of the solute plume for $R=3$, $Ra=9000$, $a/b=3$ and $f_L=0.001$. The plume is shown at dimensionless times (a) $t=0.3$, (b) $t=0.6$ and (c) $t=0.9$. These three snapshots in time show the complex nonlinear interplay between viscous and density fingering mechanisms that can lead to folding (e.g. vortex-like structure) of the plume. This observed folding dynamics is unique to the specific model configuration investigated in this work.

3. Numerical implementation

The numerical implementation of our model is similar to that presented in Bonazzi et al. (Reference Bonazzi, Morvillo, Im, Jha and de Barros2021), with the due modifications in the discretised equations to account for the concentration dependence of the solute density. To solve (2.1) numerically for cell-centred pressures, we employ a finite volume method at second-order accuracy (Ferziger, Perić & Street Reference Ferziger, Perić and Street2002). Since (2.1) is expressed in terms of the Darcy flux, i.e. the face-centred flux at the cell interface, we linearise (2.1) using the two-point flux approximation (LeVeque Reference LeVeque2002) of the Darcy flux in (2.2). In this way, we express the flux in terms of the pressure gradient, the concentration distribution and the cell interface transmissibility, which depends on viscosity. The transport problem, governed by (2.6), is solved explicitly for the concentration field, $c({\boldsymbol {x}},t)$, in terms of the advective and diffusive concentration fluxes computed at the previous time step using a second-order finite difference method.

In our numerical model we use a two-dimensional Cartesian grid where each cell has dimensions $\Delta x \times \Delta z$, with $\Delta x = 8 \times 10^{-3}$ and $\Delta z = 4 \times 10^{-3}$. The number of cells in both directions, namely $N_x$ and $N_z$, is the same, $N_x=N_z=250$. Thus, the domain lengths in the longitudinal and vertical directions are respectively $L \times H = 2 \times 1$. Our choices for $\Delta x$ and $\Delta z$ are based on a numerical grid resolution analysis that focused on quantities of interest such as the horizontally averaged concentration profile in the vertical direction and the plume's concentration variance, to ensure that the relative percentage error with respect to a more refined grid ($N_x=N_z=1000$) was below 5 %. The initial source zone in our model is an ellipse with radii $a$ and $b$ in the horizontal and vertical directions, respectively. In our paper, we investigate different ratios $a/b$ (see figure 1 and table 1) in order to explore the effects of the initial plume configuration on transport. Although the ratio $a/b$ is varied, the values of $a$ and $b$ are calculated such that the area of the ellipse, $\mathcal {S}$, at $t_0$ is the same in all the simulations; this way the same mass of solute is initially injected in the domain. For all the ratios explored, we maintain the centre of the ellipse at the dimensionless coordinate ${\boldsymbol {x}_{\boldsymbol {c}}} = (0.3,0.5)$.

Table 1.

Model input parameters. All parameter values are dimensionless following the dimensionless groups reported in Appendix A. Parameters R and Ra are varied over a range of values that are uniformly spaced between a minimum value (value on the left of the first colon) and a maximum value (value on the right of the second colon) with a constant increment or spacing between successive values (middle value between the two colons). Here $R = \ln(\mu_0/\mu_1)$ and Ra is defined in equation (2.7).

4. Comparison with previous work

In this section we compare the results originating from the numerical implementation of the governing equations as described in § 3. We test the performance of the developed numerical simulator using the same physical set-up as reported in the work of Hewitt, Neufeld & Lister (Reference Hewitt, Neufeld and Lister2013). Appendix B provides a brief description of the physical set-up and parameter values used by Hewitt et al. (Reference Hewitt, Neufeld and Lister2013). A relevant difference between our numerical simulation and that in Hewitt et al. (Reference Hewitt, Neufeld and Lister2013) is that we chose a horizontal discretisation with $\Delta x =(1024)^{-1}$, while in Hewitt et al. (Reference Hewitt, Neufeld and Lister2013) $\Delta x =(2048)^{-1}$. This was done in order to reduce the computational burden of the simulation. Another difference is that the discretisation in the vertical direction is $\Delta z =(250)^{-1}$ in Hewitt et al. (Reference Hewitt, Neufeld and Lister2013), and a vertical coordinate transformation was used to ensure high resolution of the dynamics near the fluid–fluid interface. Here, we opt for a constant vertical discretisation, $\Delta z =(750)^{-1}$. Figure 3(a) shows a comparison between the horizontally averaged concentration profiles obtained from our simulations and those reported by Hewitt et al. (Reference Hewitt, Neufeld and Lister2013) at distinct times. We can see that there is a reasonably good agreement between the two results for the considered times, although some discrepancies are present. Such discrepancies are expected given that the hydrodynamic instability, which causes the fingering behaviour, is triggered in numerical models by the introduction of a small random numerical perturbation which is amplified in time by the physics of hydrodynamic instability and the properties of the numerical method used. Differences in the amplitude and structure of the initial perturbation, which are often unreported in studies, and difference in the accuracy of the numerical discretisation method can also cause differences between studies. Figure 3(b) also illustrates a comparison of the temporal evolution of the horizontally averaged inlet solute mass flux at the top boundary of the domain. As shown in figure 3(b), our computational results are based on a single realisation of the initial perturbation, therefore displaying some noise. The mass flux estimates obtained from Hewitt et al. (Reference Hewitt, Neufeld and Lister2013) are smoother given that they represent an average from multiple realisations of the initial perturbation.

Figure 3.

(a) Comparison between the horizontally averaged concentration profiles, $\bar {c}(z,t)$, of figure 4 in Hewitt et al. (Reference Hewitt, Neufeld and Lister2013) (circles) and those obtained in our simulation (solid line), at three selected times. (b) Comparison between the mass flux $\dot {m}_s$ evolution of figure 5(b) in Hewitt et al. (Reference Hewitt, Neufeld and Lister2013) (dashed line) and that obtained in our simulation (solid line). All quantities are dimensionless in accord with the normalisation provided by Hewitt et al. (Reference Hewitt, Neufeld and Lister2013) and reported in Appendix B.

6. Effects of heterogeneity in the permeability field

Next, we focus on quantifying the effect of permeability heterogeneity on plume mixing and spreading under hydrodynamic instabilities. In order to account for heterogeneity, we consider a stratified porous medium. Stratification or layering is a common feature in many geophysical flows. For example, a layered geological formation leads to depth dependence in hydraulic properties, e.g. causing the permeability to change significantly from a clean sandstone layer of high permeability to an adjacent silty sandstone layer of low permeability. This is relevant for solute transport because it is expected that the solute plume will move faster and farther in a high-permeability layer compared with a low-permeability layer. Therefore, it is intuitive to expect enhanced plume stretching, sharper concentration gradients along the stretched interface and faster mixing across the fluid interface. We model stratification following the model described in (2.9) in Nijjer et al. (Reference Nijjer, Hewitt and Neufeld2019):

(6.1)\begin{equation} \mathrm{ln}(k) ={-}\sigma^2_{st} \cos(2{\rm \pi} nz) -\mathrm{ln}(\mathrm{I}_{0}(\sigma^2_{st})), \end{equation}

where $\sigma ^2_{st}$ is the log-permeability variance, $z$ is the vertical coordinate, ${\rm I}_0$ is the modified Bessel function of the first kind and $n$ is the wavenumber which determines the length scale of the permeability variations. In this study, we set $n=2$, thus obtaining four layers in the domain. We simulate transport for two different values of the log-permeability variance: $\sigma ^2_{st}=0.05$ and $\sigma ^2_{st}=0.1$, with the same average permeability value. We also consider the baseline case of a homogeneous medium ($\sigma ^2_{st}=0$) with a uniform permeability value equal to the average permeability in the heterogeneous cases. Regarding fluid properties, we set $R=2$ (the plume is less viscous than the ambient fluid) and $Ra=5000$ (the plume is denser than the ambient fluid). We also present results for $R=0$ for comparison. For the upcoming results, we set the inlet flux parameter to $f_L=0.01$ and consider two values of the initial source configuration: $a/b=1/3$ (vertical orientation) and $a/b=3$ (horizontal orientation).

Figure 12 shows the concentration fields for different levels of heterogeneity and values of $R$. Results are depicted at dimensionless time $t=0.1$ for $a/b=1/3$. Increasing heterogeneity in a stratified medium results in the formation of preferential flow paths in higher-permeability layers. As $\sigma ^2_{st}$ increases, the middle layer becomes more permeable, which causes additional stretching of the plume as seen in figure 12. The presence of heterogeneity changes the location of the peak concentration in the plume. For a homogeneous case ($\sigma ^2_{st}=0$), highest concentration is at the bottom of the plume because of gravitational sinking ($Ra=5000$). For layered cases, the highest concentration occurs in the high-permeability middle layer. This occurs with or without the viscosity contrast; however, viscous fingering at $R=2$ creates additional mechanisms that are superposed on the effects due to density contrast and layering. We know from the literature (Sajjadi & Azaiez Reference Sajjadi and Azaiez2013; Nicolaides et al. Reference Nicolaides, Jha, Cueto-Felgueroso and Juanes2015; Nijjer et al. Reference Nijjer, Hewitt and Neufeld2019; Bonazzi et al. Reference Bonazzi, Morvillo, Im, Jha and de Barros2021) that viscous fingering interacts with permeability heterogeneity. For example, during continuous slug injection in layered media (Nijjer et al. Reference Nijjer, Hewitt and Neufeld2019), viscous fingering dominates the mixing behaviour at early times, and when the fingers have grown to become comparable to the layer thickness scale, heterogeneity begins to dominate mixing. Such investigations of the interplay between fingering and heterogeneity mechanisms are missing for finite-size solute bodies, which we address here (see figures 12 and 13). The high flow velocity in the middle layer (figure 12e) causes tip splitting of the fingers followed by tip blunting, which we explained in our prior work; see figure 3 in Bonazzi et al. (Reference Bonazzi, Morvillo, Im, Jha and de Barros2021). The result of these interactions is rapid mixing in the high-permeability middle layer. At $R=0$, these effects are absent, and high concentrations in the middle layer persist for a longer duration.

Figure 12.

Concentration fields at $t=0.1$ for $a/b=1/3$, for three different levels of heterogeneity: (a,b) no heterogeneity ($\sigma ^2_{st}=0$), (c,d) $\sigma ^2_{st}=0.05$ and (e,f) $\sigma ^2_{st}=0.1$. (a,c,e) Results for $R=2$ and (b,d, f) the corresponding cases in the absence of viscosity contrast ($R=0$). All results are for $Ra=5000$ and $f_L=0.01$. The permeability field has been superimposed in transparency to the concentration fields; light shades correspond to higher-permeability layers, while darker shades to low-permeability ones.

Figure 13.

Maps of the velocity magnitude, $|{\boldsymbol {v}}|$, at $t=0.1$ for $f_L=0.01$ and $a/b=1/3$. The arrows indicate the velocity vector with both direction and magnitude information. The colourbars represent the logarithm of velocity magnitude, i.e. $\mathrm {log}_{10}|{\boldsymbol {v}}|$. The black line outlines the position of the plume. All results are for $Ra=5000$.

When the initial solute source is horizontally oriented ($a/b=3$, figure 14), permeability stratification does not have a significant impact on solute transport in absence of viscous fingering ($R=0$, figure 14b,d, f). This is because the source zone is not crossing multiple layers at the considered time. However, when viscosity contrast is present ($R=2$, figure 14a,c,e), stratification favours the formation of viscous fingers in the horizontal direction, a feature that was not observed for analogous cases in homogeneous media discussed in § 5. The horizontal finger formation can be explained by observing figure 15(a,c,e). As the difference in velocity magnitude between layers increases, the upper and lower edges of the plume experience lower velocity magnitudes compared with the centre. The horizontal direction of the velocity field causes the centre of the plume to travel faster than its edges and split into two fingers, following the tip splitting mechanism. At $\sigma ^2_{st}=0.05$ (figure 15c), the upper finger eventually merges into the dominating lower finger, following the finger coarsening mechanism (Jha et al. Reference Jha, Cueto-Felgueroso and Juanes2011b). At $\sigma ^2_{st}=0.1$, heterogeneity is strong enough to sustain the growth of both fingers for a longer time (figure 14e). Eventually, the fingers will merge as in the case of a vertically oriented source ($a/b=1/3$), and permeability heterogeneity will be the dominant agent affecting solute mixing and transport.

Figure 14.

Concentration fields for $a/b=3$ at $t=0.1$ for three levels of heterogeneity: (a,b) $\sigma ^2_{st}=0$ (no heterogeneity), (c,d) $\sigma ^2_{st}=0.05$ and (e, f) $\sigma ^2_{st}=0.1$. All results are for $Ra=5000$ and $f_L=0.01$.

Figure 15.

Maps of the velocity magnitude, $|{\boldsymbol {v}}|$, for $a/b=3$ at $t=0.1$ for $f_L=0.01$. The arrows indicate the velocity vector with both direction and magnitude information. The colourbars represent the logarithm of velocity magnitude, i.e. $\mathrm {log}_{10}|{\boldsymbol {v}}|$. The black line outlines the position of the plume. All results are for $Ra=5000$.

In this section, the layer thickness is comparable to the size of the initial source. The horizontally oriented source is inside the higher-permeability layer, i.e. the permeability values in that layer (albeit not constant due to the sinusoidal nature of the field) are higher compared with the mean permeability. However, we can see that even such slight variation of the permeability is able to affect the transport behaviour. This suggests that the structure of the permeability field might impact the interplay between density and viscosity contrast in driving mixing. If the permeability field structure is reversed (i.e. switching the locations of high-permeability layers with the locations of the low-permeability layers), or if the solute injection is located between layers, the impact of the structure of the permeability field on mixing might be different.

7. Summary

The need for modelling the effects of viscosity and density contrasts in the spreading and mixing dynamics of a solute body has been recognised and addressed in the literature given its importance in contaminant transport, energy recovery and storage in the subsurface environment. In all these applications, solutes are characterised by viscosity and density values that differ from those of the ambient fluid. Furthermore, these solutes are injected into the subsurface through source zones of different spatial configurations. Our work focuses on systematically investigating the impact of the initial configuration of the source zone on the transport behaviour in the presence of hydrodynamic instabilities. The numerical analysis carried out in the present work illustrates the importance of both the initial solute source configuration and the intensity of the inlet (background) flux in controlling the relative significance of $R$ and $Ra$ in the mixing and spreading behaviour. In the following, we highlight the key findings of our analysis.

Effect of viscosity contrast on mixing: For a homogeneous medium and low inlet or background flux, gravitational fingering is impeded when the solute is more viscous than the ambient fluid, i.e. $R<0$. When $R>0$, viscous fingering occurs in both horizontal and vertical directions and in the latter case, viscous fingering and gravitational fingering effects overlap. For $R<0$, better mixing is achieved with a vertically oriented initial source. If more mixing is desired under $R>0$ conditions, a horizontally oriented solute source is preferable at a low background flux, whereas a vertical source is preferable at high background flux. At higher inlet flux, when $R>0$, a horizontally oriented source does not result in the best degree of mixing because the dominance of viscous fingering (as the primary driver of mixing) allows a vertically extended source to form a larger number of fingers. The importance of viscosity contrast in enhancing mixing when the inlet flux is higher means that the degree of mixing is lower when no viscosity contrast is present ($R\approx 0$).

Effect of density contrast on mixing: For a homogeneous medium and low inlet flux, higher values of $Ra$ were also found to allow lower rates of diffusive mass transfer. The impact of $Ra$ on mixing was found to be almost negligible when $R>0$. A horizontally oriented initial source zone allows gravitational fingers to form and grow in the vertical direction. The number of fingers increases with $Ra$ in accordance with the linear stability analysis theory. At higher inlet flux, the role of $Ra$ in mixing is minor for a vertically oriented source for any value of $R$ considered. For a horizontally oriented source, $Ra$ appears to affect mixing when $R<0$.

Effect of background flow velocity on mixing: We established that a higher mean flow velocity hampers gravitational fingering, thus rendering viscous fingering the dominant force driving mixing. As a consequence, the source orientation becomes a significant factor when $R<0$. Under these conditions, we observed the following: (1) for $a/b=1/3$ (vertically oriented source), a higher flow velocity allows the formation of viscous fingers, which will lead to enhanced mixing; while (2) for $a/b=3$ (horizontally oriented source) the absence of viscous fingering and the low velocity inside the plume result in poor mixing.

Relationship between mixing and plume shape: The numerical analysis presented here allowed us to investigate the relationship between the degree of mixing and the length of the interface between the solute and the ambient fluid (i.e. plume shape). Traditionally, one would expect enhanced mixing rates for increased elongation of the plume's interface with the surrounding fluid. In the context of aquifer remediation, hydrologists and engineers aim to maximise the interface between the contaminant and groundwater in order to enhance dilution and reaction rates (e.g. Bagtzoglou & Oates Reference Bagtzoglou and Oates2007; Kitanidis & McCarty Reference Kitanidis and McCarty2012; Trefry et al. Reference Trefry, Lester, Metcalfe, Ord and Regenauer-Lieb2012; Piscopo, Neupauer & Mays Reference Piscopo, Neupauer and Mays2013). As shown in § 5.2.1, the relationship between mixing and plume shape is non-trivial when fingering instabilities are present. The results in § 5.2.1 shed novel insights on the role of source orientation, inlet flux and fingering instabilities (i.e. viscosity and density contrasts) in controlling the plume shape and its mixing dynamics. Our results show that the same degree of mixing can be achieved for different plume shapes.

Effect of permeability heterogeneity on mixing: To account for the heterogeneity, we considered a stratified permeability field. We found that, for a vertically oriented solute source, heterogeneity in the porous medium affects the location of the plume's highest concentration values by overlaying the effect of gravitational instability and favouring higher concentration values in the higher-permeability layer. The high flow velocity in the more permeable layer causes tip splitting and results in rapid mixing. For a horizontally oriented solute source, the combined effect of viscosity contrast and permeability heterogeneity results in fingering along the flow direction, a result not observed for the homogeneous media analysis presented in § 5.

The simulations carried out in this work highlight the importance of the initial solute source configuration and the intensity of the inlet flux in modulating both the plume's spreading and mixing behaviour in homogeneous and heterogeneous permeability fields. Additional research is needed to understand these complex interactions in three-dimensional settings, larger computational domains, in the absence of interaction with the bottom boundary of the domain, with different injection modes and distinct conceptualisations of the permeability field, including higher levels of permeability heterogeneity.