2.1 Non-dimensionalization

We non-dimensionalize the governing equations by the width of the porous medium $a$ , the injection flux $Q$ , the mean permeability $\bar{k}$ , viscosity $\unicode[STIX]{x1D707}_{2}$ , pressure $\unicode[STIX]{x1D707}_{2}Q/\bar{k}$ and convective time scale $\unicode[STIX]{x1D719}a^{2}/Q$ . We also change to a reference frame moving with the average speed of the injected fluid by setting $\tilde{x}=x-Qt/a$ and $\tilde{u} =u-Q/a$ . Equations (2.1)–(2.4) become

(2.5a,b ) $$\begin{eqnarray}-(\tilde{u} +1)\unicode[STIX]{x1D707}=k\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\tilde{x}},\quad -v\unicode[STIX]{x1D707}=k\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y},\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{u} }{\unicode[STIX]{x2202}\tilde{x}}+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}=0, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}+\tilde{u} \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\tilde{x}}+v\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}=\frac{1}{Pe}\left(\frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}\tilde{x}^{2}}+\frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}y^{2}}\right), & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}(c)=\text{e}^{-Rc}. & \displaystyle\end{eqnarray}$$

The system is described by two non-dimensional parameters, the Péclet number $Pe=Q/D$ , which characterizes the strength of advection to diffusion and the log-viscosity ratio $R$ , as well as the non-dimensional permeability, $k(y)$ . For notational convenience we drop the tildes in all subsequent expressions.

We consider three different cases for the log-viscosity ratio: the neutrally stable case $R=0$ , where the injected and ambient fluids have the same viscosity; the ‘stable’ case $R<0$ , where the injected fluid is more viscous than the ambient; and the ‘unstable’ case $R>0$ , where the injected fluid is less viscous than the ambient. The Péclet number provides a measures of the relative strength of advection and diffusion, and in this work we will focus on the limit of large, but finite, Péclet number. This is the typical limit in most geologic scenarios.

2.2 Choice of permeability

In this paper, we consider only layered heterogeneous media for which $k=k(y)$ . In fact, for all the numerical results presented here, we further restrict our attention to log-permeabilities which vary sinusoidally (De Wit & Homsy Reference De Wit and Homsy1997a ,Reference De Wit and Homsy b ; Sajjadi & Azaiez Reference Sajjadi and Azaiez2013), such that

(2.9) $$\begin{eqnarray}\ln (k)=-\unicode[STIX]{x1D70E}\cos (2\unicode[STIX]{x03C0}ny)-\text{ln}\,(\text{I}_{0}(\unicode[STIX]{x1D70E})),\end{eqnarray}$$

where $\text{I}_{0}$ is the modified Bessel function of the first kind, which ensures a unit average dimensionless permeability. This simplification retains the dominant physics of permeability heterogeneities in the form of layering, while being able to be described by two parameters instead of the infinite space of possible permeability functions. These two parameters are the log-permeability variance $\unicode[STIX]{x1D70E}$ and wavenumber $n$ , which measure the strength and inverse of the length scale of the permeability variations, respectively. While some of our results are presented for general $k(y)$ , all of the numerical simulations in this paper use this form for the permeability structure.

In the absence of hydrodynamic instabilities, as described in §§ 3–5, there is no mechanism for dynamic interactions between layers and so $n$ can be scaled out of the system. This is done by introducing rescaled variables ${\hat{y}}=ny$ , $\hat{x}=nx$ , $\hat{t}=n^{2}t$ , in which case the flow evolves exactly as it would with $n=1$ , but with an effective Péclet number $\hat{Pe}=Pe/n^{2}$ . For clarity, we therefore limit our analysis in §§ 3–5 to the case where $n=1$ . If however the flow is unstable, as in § 6, there can be a competition between the evolving wavelength of the viscous fingering and the imposed wavelength of the permeability structure, resulting in rich intermediate-time dynamics for which the value of $n$ can be important. We therefore explicitly consider the dependence of the flow on the number of layers in § 6.

2.3 Boundary and initial conditions

Following previous work (e.g. Tan & Homsy Reference Tan and Homsy1986; Nijjer et al. Reference Nijjer, Hewitt and Neufeld2018), we consider flows that are periodic in the transverse ( $y$ ) direction,

(2.10) $$\begin{eqnarray}c(x,0,t)=c(x,1,t),\end{eqnarray}$$

(2.11a,b ) $$\begin{eqnarray}u(x,0,t)=u(x,1,t),\quad v(x,0,t)=v(x,1,t).\end{eqnarray}$$

The upstream and downstream fluxes are zero (in the moving frame) and the transverse velocity vanishes in the far field so that,

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{1}u\,\text{d}y\rightarrow 0\quad \text{as }x\rightarrow \pm \infty , & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}\rightarrow 0\quad \text{as }x\rightarrow \pm \infty , & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle v\rightarrow 0\quad \text{as }x\rightarrow \pm \infty . & \displaystyle\end{eqnarray}$$

We initialize the concentration field to have a step jump,

(2.15) $$\begin{eqnarray}c(x,t=0)=c_{0}(x)=H(-x),\end{eqnarray}$$

where $H(x)$ is the Heaviside function.

2.4 Diagnostic quantities

In order to investigate how the macroscopic features of the flow evolve we focus, in the following analysis, on the evolution of the transversely averaged concentration $\overline{c}(x,t)=\int _{0}^{1}c\,\text{d}y$ and the mixing length $h(t)$ , defined below. To understand how $\overline{c}$ evolves, we start by decomposing the advection–diffusion equation (2.7), into two coupled equations for the transversely averaged concentration and the deviations $c^{\prime }=c-\overline{c}$ ,

(2.16) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\overline{uc^{\prime }}}{\unicode[STIX]{x2202}x}=\frac{1}{Pe}\frac{\unicode[STIX]{x2202}^{2}\overline{c}}{\unicode[STIX]{x2202}x^{2}},\end{eqnarray}$$

and

(2.17) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c^{\prime }}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}uc^{\prime }}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}u\overline{c}}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}\overline{uc^{\prime }}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}vc^{\prime }}{\unicode[STIX]{x2202}y}=\frac{1}{Pe}\left(\frac{\unicode[STIX]{x2202}^{2}c^{\prime }}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}c^{\prime }}{\unicode[STIX]{x2202}x^{2}}\right),\end{eqnarray}$$

where $\overline{f}=\int _{0}^{1}f\,\text{d}y$ denotes the transverse average of the quantity $f$ . These equations are obtained by averaging the advection–diffusion equation in the transverse direction and by subtracting the averaged equation from the unaveraged one. We return to this decomposition when describing the evolution of the transversely averaged concentration.

The region over which the two fluids have spread, the ‘mixing length’ $h(t)$ , is a measure of the streamwise extent of interpenetration of the two fluids. The mixing length provides a global measure of dispersion or spreading, and, while other quantities can give more direct measurements of the total amount of mixing (e.g. the scalar dissipation rate; Jha et al. Reference Jha, Cueto-Felgueroso and Juanes2011a ), the mixing length provides a clear and physically useful measure of the region over which the two fluids have spread and are mixing. Here we define the mixing length to be the variance in concentration about the initial condition $c_{0}(x)$ (2.15),

(2.18) $$\begin{eqnarray}h(t)=\sqrt{\displaystyle \frac{\displaystyle \int _{-\infty }^{\infty }x^{2}(\overline{c}-c_{0})^{2}\,\text{d}x}{\displaystyle \int _{-\infty }^{\infty }(\overline{c}-c_{0})^{2}\,\text{d}x}}.\end{eqnarray}$$

Note that we use this measure instead of the more common definition, $h^{\ast }=x|_{\overline{c}=0.01}-x|_{0.99}$ , because it more accurately captures the spreading behaviour when the concentration field has long tails (see appendix A).

Throughout this paper we make reference to the interface between the two fluids. Since the two fluids are fully miscible, there is no precise boundary. Instead, where we refer to the interface, we loosely mean the region around the $c=1/2$ contour over which the concentration varies significantly.

2.5 Numerical method

We briefly summarize the numerical method here, but avoid a detailed exposition as the approach is very similar to that used by Nijjer et al. (Reference Nijjer, Hewitt and Neufeld2018). We start by combining (2.5), (2.6) and (2.8) and writing the velocity in terms of a streamfunction, $(u,v)=(\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}/\unicode[STIX]{x2202}y,-\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}/\unicode[STIX]{x2202}x)$ , to give

(2.19) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}y^{2}}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}x}\left(R\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}\right)-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}y}\left(R\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}+\frac{\text{d}\ln k}{\text{d}y}\right)=R\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}+\frac{\text{d}\ln k}{\text{d}y}.\end{eqnarray}$$

The transverse velocity boundary conditions (2.11) become

(2.20) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}(x,0,t)=\unicode[STIX]{x1D6F9}(x,1,t).\end{eqnarray}$$

We imposed the boundary conditions in the flow direction (2.12)–(2.14), at a finite distance $x=\pm \unicode[STIX]{x1D6E4}$ so that,

(2.21) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}x}=0\quad \text{at }x=\pm \unicode[STIX]{x1D6E4},\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}$ was chosen to be sufficiently large such that the mixing region was always far from the boundaries. In fact, we increased $\unicode[STIX]{x1D6E4}$ over the course of each simulation so that it grew as the mixing length grew, ensuring that the mixing region was always far from the boundaries while the fine-scale dynamics at early times remained resolved. The domain was discretized on an adaptive rectangular grid which coarsened over time as the concentration gradients weakened.

The concentration field was initialized to an almost sharp interface, to avoid any numerical instabilities, by setting

(2.22) $$\begin{eqnarray}c_{0}={\textstyle \frac{1}{2}}-{\textstyle \frac{1}{2}}\text{erf}\left(x/\sqrt{t_{0}}\right)+r(x,y)\text{e}^{-x^{2}/t_{0}}.\end{eqnarray}$$

Here $t_{0}$ is a small time origin, which we set to $t_{0}=10^{-6}$ in all simulations, and $r$ corresponds to a uniformly distributed random function, on the interval $[0,10^{-4}]$ , which was added to help trigger any instabilities.

At each time step (2.19) was solved using a multi-grid solver Adams (Reference Adams1999). Equation (2.3) was advanced in time using a third-order Runge–Kutta scheme. All spatial derivatives were discretized using sixth-order compact finite differences except for the advection term in (2.3), $\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c$ , which was discretized using a third-order upwinding scheme. We validated our scheme by comparison to the linear stability analysis of Pramanik & Mishra (Reference Pramanik and Mishra2015a ), and we chose spatial discretizations such that the smallest scales of the flow were always well resolved. In particular, we ensured that in cases when the interface was unstable to small-scale fingering the resolution was sufficiently high that doubling the grid spacing had no qualitative effect on the dynamics of the flow or diagnostic quantities.

Figure 2. Evolution of the concentration field for $R=0$ and $(\unicode[STIX]{x1D70E},Pe,n)=(1,100,1)$ . (a) The imposed permeability $k(y)=\text{e}^{-\text{cos}(2\unicode[STIX]{x03C0}y)}/\text{I}_{0}(1)$ . (b) Evolution of the mixing length, $h$ , as a function of time, $t$ . The dots correspond to the snapshots (c–e). (c–e) Plots of the concentration field with overlain streamlines versus $x/h$ and $y$ at (c) $t=5\times 10^{-4}$ , (d) $t=0.32$ and (e) $t=200$ .

3.1 Early-time behaviour: initial diffusion

At early times, the streamwise concentration gradient between the fluids is large and the concentration is transversely homogeneous. In this case, diffusion across the interface dominates and the primary balance in the advection–diffusion equation is

(3.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}t}=\frac{1}{Pe}\frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x^{2}}.\end{eqnarray}$$

Using the initial and boundary conditions in § 2.3 the concentration evolves self-similarly as

(3.4) $$\begin{eqnarray}c=\overline{c}=\frac{1}{2}+\frac{1}{2}\text{erf}\left(-\frac{x}{\sqrt{4t/Pe}}\right),\end{eqnarray}$$

which holds at all times when the permeability is homogeneous ( $k=1$ ). The mixing length grows like $h\sim t^{1/2}$ and can be calculated explicitly by substituting (3.4) into (2.18) (figure 3 a).

This behaviour always holds initially, irrespective of the parameter choices, since diffusive growth of the interface $O(t^{1/2}/Pe^{1/2})$ always outpaces advective spreading $O(\unicode[STIX]{x0394}ut)$ (where $\unicode[STIX]{x0394}u$ is a characteristic spreading velocity). In fact, the transition to the intermediate regime occurs precisely when the growth rates become equal, giving a transition time $t=O(1/Pe\unicode[STIX]{x0394}u^{2})$ .

3.2 Intermediate-time behaviour: advection

After a time $O(1/Pe\unicode[STIX]{x0394}u^{2})$ , spreading induced by the difference in permeability overtakes longitudinal diffusion. The leading-order balance in (2.3) becomes

(3.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}=0,\end{eqnarray}$$

and the flow simply stretches the diffused solution that arises from the early-time regime. In fact, since the rate of advective stretching is much faster than diffusion to good approximation, we can ignore the effects of the early-time regime completely and the solution to (3.5) is simply the travelling wave

(3.6) $$\begin{eqnarray}c=c_{0}(x-u(y)t)=H(u(y)-x/t),\end{eqnarray}$$

given the initial condition (2.15). The transversely averaged concentration, $\overline{c}(x,t)$ can be calculated by averaging (3.6),

(3.7) $$\begin{eqnarray}\overline{c}(x,t)=\int _{0}^{1}H(u(y)-x/t)\,\text{d}y.\end{eqnarray}$$

The model gives good agreement with the numerical simulations (figure 3 b) and is able to reproduce the asymmetric profiles, which arise due to the fact that the permeability is not symmetric about $k=1$ and collapses upon rescaling the streamwise coordinate by $t$ . Since the interface is stretched at a constant rate, the mixing length grows like $h\sim \unicode[STIX]{x0394}ut$ (figures 3 a and 4 b).

Figure 3. Evolution of the mixing length and transversely averaged concentration for $R=0$ in the intermediate-time regime. (a) Scaled plot of the mixing length $h(t)$ for $(Pe,\unicode[STIX]{x1D70E})=(500,1)$ , $n=\{1,2,3,4,6,8\}$ and $(Pe,n)=(500,1)$ , $\unicode[STIX]{x1D70E}=\{0.1,0.2,0.4,0.6,0.8,1,1.2,1.4,1.8,2\}$ and $(\unicode[STIX]{x1D70E},n)=(1,1)$ , $Pe=\{1,2,3,4,5,6,8,15,20\}\times 10^{2}$ . The black lines correspond to the predictions for the mixing length calculated from the analytical solutions (3.4) (dashed) and (3.7) (dotted). (b) Plot of the transversely averaged concentration, $\overline{c}(x,t)$ versus $x/t$ for $(R,Pe,n)=(0,500,1)$ , $\unicode[STIX]{x1D70E}$ ranging from $0.2$ to $1.8$ and ten logarithmically spaced times between $1$ and $3$ . The characteristic spreading velocity is taken to be the maximum velocity difference between the layers, $\unicode[STIX]{x0394}u\sim \sinh (\unicode[STIX]{x1D70E})/\text{I}_{0}(\unicode[STIX]{x1D70E})$ .

3.3 Late-time behaviour: shear-enhanced dispersion

The transition to the late-time regime occurs once the concentration has diffused across the entire channel, homogenizing the concentration in the transverse direction; this occurs at a time $O(Pe)$ . In this case, the mixing zone is long and thin and transverse diffusion balances longitudinal advection (cf. Taylor dispersion e.g. Taylor Reference Taylor1953; Aris Reference Aris1956). This is in contrast to the previous regime, when the flow evolved purely by longitudinal advection. In the limit of small deviations from the mean $c^{\prime }\ll \overline{c}$ and a long, thin mixing zone, (2.17) reduces to

(3.8) $$\begin{eqnarray}(k-1)\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}x}=\frac{1}{Pe}\frac{\unicode[STIX]{x2202}^{2}c^{\prime }}{\unicode[STIX]{x2202}y^{2}},\end{eqnarray}$$

while the transversely averaged concentration still evolves according to (2.16). Given that $\overline{c}$ is independent of $y$ , we integrate this equation twice and impose periodicity and zero-mean deviations ( $\int _{0}^{1}c^{\prime }=0$ ), to give

(3.9) $$\begin{eqnarray}c^{\prime }=Pe\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}x}\left[\int _{0}^{y}\int _{0}^{\unicode[STIX]{x1D701}}(k(\unicode[STIX]{x1D702})-1)\,\text{d}\unicode[STIX]{x1D702}\,\text{d}\unicode[STIX]{x1D701}-\int _{0}^{1}\int _{0}^{s}\int _{0}^{\unicode[STIX]{x1D701}}(k(\unicode[STIX]{x1D702})-1)\,\text{d}\unicode[STIX]{x1D702}\,\text{d}\unicode[STIX]{x1D701}\,\text{d}s\right].\end{eqnarray}$$

Substituting and solving for the convective flux in (2.16), using the expression for the velocity (3.1), leads to

(3.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\overline{uc^{\prime }}}{\unicode[STIX]{x2202}x}=-Pe\frac{\unicode[STIX]{x2202}^{2}\overline{c}}{\unicode[STIX]{x2202}x^{2}}\int _{0}^{1}\left[\int _{0}^{y}(k(\unicode[STIX]{x1D702})-1)\,\text{d}\unicode[STIX]{x1D702}\right]^{2}\,\text{d}y.\end{eqnarray}$$

This convective flux can be written in the form of an effective diffusivity such that (2.16) reduces to

(3.11a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\frac{1}{Pe^{\ast }}\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}x}\right],\quad \frac{1}{Pe^{\ast }}=\frac{1}{Pe}(1+Pe^{2}{\mathcal{S}});\end{eqnarray}$$

where

(3.12) $$\begin{eqnarray}{\mathcal{S}}=-\frac{1}{Pe\,\unicode[STIX]{x2202}\overline{c}/\unicode[STIX]{x2202}x}\int _{0}^{1}uc^{\prime }\,\text{d}y=\int _{0}^{1}\left[\int _{0}^{y}(k(\unicode[STIX]{x1D702})-1)\,\text{d}\unicode[STIX]{x1D702}\right]^{2}\,\text{d}y\end{eqnarray}$$

is the shear-enhanced dispersivity, which only depends on the permeability structure (cf. Van den Broeck & Mazo Reference Van den Broeck and Mazo1983).

Figure 4. Evolution of the mixing length and transversely averaged concentration for $R=0$ in the late-time regime. (a) Shear-enhanced dispersivity ${\mathcal{S}}$ versus $\unicode[STIX]{x1D70E}$ with asymptotic limits (B 3) (dashed) and (B 6) (dot-dashed). (b) Scaled plot of $h$ versus $t$ for $(Pe,\unicode[STIX]{x1D70E})=(500,1)$ , $n=\{1,2,3,4,6,8\}$ and $(Pe,n)=(500,1)$ , $\unicode[STIX]{x1D70E}=\{0.1,0.2,0.4,0.6,0.8,1,1.2,1.4,1.8,2\}$ and $(\unicode[STIX]{x1D70E},n)=(1,1)$ , $Pe=\{1,2,3,4,5,6,8,15,20\}\times 10^{2}$ . The black lines correspond to the predictions for the mixing length calculated from the analytical solutions (3.7) (dotted) and (3.4) with (3.11) (dashed). (c) Plot of the transversely averaged concentration, $\overline{c}(x,t)$ versus the late-time similarity variable $x/\sqrt{t/Pe^{\ast }}$ for the same parameters as (b) at $t=200$ . The theoretical solution (3.4) with (3.11) is given by the dashed line.

For our choice of sinusoidally varying log-permeability, this integral cannot be solved analytically, but is instead integrated numerically for varying $\unicode[STIX]{x1D70E}$ and plotted in figure 4(a). In appendix B, we derive the asymptotic limits of the shear-enhanced dispersivity for large and small $\unicode[STIX]{x1D70E}$ (given as dashed and dot-dashed lines respectively in figure 4 a).

The solution to (3.11a ) is again the similarity solution (3.4), but now with a modified Péclet number $Pe^{\ast }$ (3.11b ) (see figure 4 c). When the total dispersivity is dominated by the shear-enhanced dispersivity, $Pe^{2}{\mathcal{S}}\gg 1$ , the effective dispersion scales like $Pe^{\ast }\sim 1/(Pe\unicode[STIX]{x0394}u^{2})$ . In figure 4(b) we use this scaling to collapse the mixing length as a function of time over a range of parameters.

In summary, in the presence of permeability layering but in the absence of viscosity variations, the flow evolves through three regimes: early-time diffusion, intermediate-time advection and late-time shear-enhanced dispersion.

4.1 Vertical flow equilibrium

Unlike when the viscosities are equal, $R=0$ , we cannot simply integrate (2.5) to give a fixed expression for the velocity. As noted earlier, we find that when the injected fluid is more viscous than the ambient, $R<0$ , the streamwise velocity is reduced, whereas when the injected fluid is less viscous than the ambient, $R>0$ , the streamwise velocity is increased. Under the assumption that the flow is long and thin, the pressure is only a function of the longitudinal coordinate and is constant along any transverse slice to leading order. This limit, often referred to as ‘vertical flow equilibrium’ (Yortsos Reference Yortsos1995), implies that

(4.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D707}(u+1)}{k}=-\frac{\text{d}p}{\text{d}x}\approx \text{const}.\end{eqnarray}$$

Combining (4.1) with the fact that the flux vanishes in any transverse slice in the moving frame, the velocity can be written as

(4.2) $$\begin{eqnarray}u(x,y,t)=\frac{k(y)}{\unicode[STIX]{x1D707}(x,y,t)}\left[\int _{0}^{1}\frac{k(s)}{\unicode[STIX]{x1D707}(x,s,t)}\,\text{d}s\right]^{-1}-1.\end{eqnarray}$$

If the viscosity is uniform, then the permeability sets the velocity, $u=k-1$ , as in (3.1). If, instead, the permeability is uniform, then the viscosity sets the velocity, $u(y)=\unicode[STIX]{x1D707}^{-1}/\int _{0}^{1}\unicode[STIX]{x1D707}^{-1}\,\text{d}y-1$ (this leads to the fast low-viscosity fingers and slow high-viscosity fingers characteristic of the viscous-fingering instability). When both the permeability and viscosity vary, depending on the sign of $R$ and $c^{\prime }$ , the permeability and viscosity can interact either constructively or destructively. The effect of varying viscosity is only important at the interface; far upstream and downstream, where the viscosity is uniform, the velocity variations are simply imposed by the structure of the permeability field.

Decomposing the concentration into the transverse average and deviations $c=\overline{c}(x)+c^{\prime }(x,y)$ , equation (4.2) becomes independent of the average concentration and only depends on the transverse variations,

(4.3) $$\begin{eqnarray}u(x,y,t)=\frac{k(y)}{\text{e}^{-Rc^{\prime }(x,y,t)}}\left[\int _{0}^{1}\frac{k(s)}{\text{e}^{-Rc^{\prime }(x,s,t)}}\,\text{d}s\right]^{-1}-1.\end{eqnarray}$$

In the case of sinusoidal log-permeability variations, equation (4.3) is

(4.4) $$\begin{eqnarray}u(x,y,t)=\frac{\text{e}^{-\unicode[STIX]{x1D70E}\cos (2\unicode[STIX]{x03C0}y)+Rc^{\prime }(x,y,t)}}{\displaystyle \int _{0}^{1}\text{e}^{-\unicode[STIX]{x1D70E}\cos (2\unicode[STIX]{x03C0}s)+Rc^{\prime }(x,s,t)}\,\text{d}s}-1.\end{eqnarray}$$

4.2 Intermediate-time behaviour: viscously coupled advection

After the early-time diffusion regime the flow transitions to the intermediate-time regime dominated by advective spreading. The effect of a non-zero viscosity ratio on the evolution of the mixing length is shown in figure 6(a). Similar to the $R=0$ case, the mixing zone grows linearly in time, however as $R$ is increased, the growth rate of the mixing zone, ${\dot{h}}=\text{d}h/\text{d}t$ , increases. The nearly uniform spacing between the curves suggests that the growth rate varies linearly in $R$ . The transversely averaged concentration is again asymmetric and evolves self-similarly (figure 6 b), although it differs appreciably from the $R=0$ case at the downstream tips (cf. snapshots in figure 5 a,c,e).

Figure 6. Evolution of the mixing length and transversely averaged concentration for $|R|<\unicode[STIX]{x1D70E}$ in the intermediate-time regime. (a) Evolution of the mixing length for $(\unicode[STIX]{x1D70E},Pe,n)=(1,2000,1)$ and $R$ ranging from $-0.4$ to $0.4$ ; the inset shows the same data normalized by the neutral displacement mixing length. (b) Plot of $\overline{c}$ versus $x/t$ for the same parameters as (a) for $10$ logarithmically spaced times in the range $1\leqslant t\leqslant 3$ . (c) Plot of the concentration deviations $c^{\prime }$ as a function $y$ at three different points in $x$ corresponding to $\overline{c}=0.75$ (red), $\overline{c}=0.5$ (green) and $\overline{c}=0.25$ (blue) for $(R,\unicode[STIX]{x1D70E},Pe,n)=(-0.4,1,2000,1)$ . The longitudinally averaged concentration deviations (averaged over the mixing zone), $\int _{h}c^{\prime }\text{d}x/h$ , is given by the dashed black line. (d) Plot of the spreading rate ${\dot{h}}$ calculated by least-squares fitting a function of the form $h=h_{0}+{\dot{h}}t$ to the numerical results for $t$ in the range $1\leqslant t\leqslant 3$ , for $Pe=2000$ (circles) and $Pe=4000$ (squares). The theoretical predictions for ${\dot{h}}$ , calculated using (4.5) and (3.7), are given by the solid lines.

In this regime, we first note that diffusion is negligible. This results in concentration deviations that are almost exactly either $c^{\prime }=-\overline{c}$ or $c^{\prime }=1-\overline{c}$ (figure 6 c). To estimate the overall effect of the deviations on the velocity, we average the deviations across the length of the fingered region, which leads to a roughly sinusoidal variation across the domain aligned with the permeability structure and with magnitude $\simeq 1/2$ (dashed black line in figure 6 c). Substituting these average deviations into (4.4), the mean streamwise velocity reduces to

(4.5) $$\begin{eqnarray}u=\frac{\text{e}^{-(\unicode[STIX]{x1D70E}+R/2)\cos (2\unicode[STIX]{x03C0}y)}}{\displaystyle \int _{0}^{1}\text{e}^{-(\unicode[STIX]{x1D70E}+R/2)\cos (2\unicode[STIX]{x03C0}s)}\,\text{d}s}-1.\end{eqnarray}$$

This approximate model results in intermediate-time dynamics equivalent to the neutrally stable case, but with an effective log-permeability ratio $\unicode[STIX]{x1D70E}_{eff}=\unicode[STIX]{x1D70E}+R/2$ . Given (4.5), we can calculate $h$ as in § 3.2, and extract ${\dot{h}}$ by fitting a linear profile $h={\dot{h}}t$ . Modelling the effective permeability in this way gives reasonably good agreement with the numerical simulations (figure 6 d), although it underestimates the spreading rate for large $|R|$ . This is because, at the boundary of the forward-propagating and backward-propagating tips, the velocity is faster, $u=\text{e}^{\unicode[STIX]{x1D70E}+R}/\text{I}_{0}(\unicode[STIX]{x1D70E})$ , and slower, $u=\text{e}^{-(\unicode[STIX]{x1D70E}+R)}/\text{I}_{0}(\unicode[STIX]{x1D70E})$ , respectively, than this model predicts.

4.3 Late-time behaviour: viscosity-dependent shear-enhanced dispersion

At late times, after advectively spreading, the concentration evolves diffusively again. This evolution is analogous to the late-time behaviour of the neutrally stable case but the addition of viscosity variations modifies the effective diffusivity.

As in § 3.3, we assume that the flow is long and thin, transverse velocities are negligible so the fluid flow is predominantly in the streamwise direction and the concentration deviations are small and evolve on a much faster time scale than their transverse average. Equation (2.16) remains unchanged, but (3.8) becomes

(4.6) $$\begin{eqnarray}\left(\frac{k\text{e}^{Rc^{\prime }(x,y,t)}}{\displaystyle \int _{0}^{1}k\text{e}^{Rc^{\prime }(x,s,t)}\,\text{d}s}-1\right)\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}x}Pe=\frac{\unicode[STIX]{x2202}^{2}c^{\prime }}{\unicode[STIX]{x2202}y^{2}}(x,y,t)\end{eqnarray}$$

because of the dependence of the velocity on the concentration through (4.3). Again we impose periodic boundary conditions and zero-mean deviations.

Before solving, we first rescale the deviations, $\tilde{c}=Rc^{\prime }$ , such that (4.6) becomes

(4.7) $$\begin{eqnarray}\left(\frac{k\text{e}^{\tilde{c}}}{\displaystyle \int _{0}^{1}k\text{e}^{\tilde{c}}\,\text{d}s}-1\right)\unicode[STIX]{x1D6FD}=\frac{\unicode[STIX]{x2202}^{2}\tilde{c}}{\unicode[STIX]{x2202}y^{2}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}=R\,Pe\,\unicode[STIX]{x2202}\overline{c}/\unicode[STIX]{x2202}x=Pe\,\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D707}(\overline{c})/\unicode[STIX]{x2202}x$ , incorporates all of the parameters in the problem, and can be thought of as a rescaled bulk concentration or viscosity gradient. Since $\unicode[STIX]{x1D6FD}$ is independent of $y$ , we can solve for $\tilde{c}$ by integrating (4.7) twice. We perform this integration numerically, although the limits of small $\unicode[STIX]{x1D6FD}$ , which is relevant here and is considered in § 4.3.1, and large $\unicode[STIX]{x1D6FD}$ , which we find to be useful and is considered later in § 5.1, can be treated analytically. Having solved for $\tilde{c}$ , we then calculate the shear-enhanced dispersivity,

(4.8) $$\begin{eqnarray}{\mathcal{S}}=-\frac{1}{\unicode[STIX]{x1D6FD}}\int _{0}^{1}u\tilde{c}\,\text{d}y,\end{eqnarray}$$

(cf. (3.12)). Solutions of ${\mathcal{S}}$ for $\unicode[STIX]{x1D70E}=1$ and $R>0$ and $R<0$ are given in figure 7. When $R=0$ , (4.7) reduces to (3.8) and ${\mathcal{S}}$ is exactly as described by (3.12), a permeability-dependent constant. When the viscosity varies, the dispersivity is enhanced when the injected fluid is less viscous ( $R>0$ ) and diminished when the injected fluid is more viscous ( $R<0$ ), than the ambient fluid. Since $\unicode[STIX]{x1D6FD}\propto \unicode[STIX]{x2202}\overline{c}/\unicode[STIX]{x2202}x$ , and diffusion always causes concentration gradients to diminish in time, $\unicode[STIX]{x1D6FD}$ generally decreases over time. When $|\unicode[STIX]{x1D6FD}|$ is large (i.e. the viscosity gradient is large, as at early times), and $R>0$ , ${\mathcal{S}}$ diverges, whereas when $R<0$ , ${\mathcal{S}}$ tends to zero. The former limit is unphysical and corresponds to scenarios where the interface is unstable. The latter case is discussed in more detail in § 5, in the context of stable injections. When $|\unicode[STIX]{x1D6FD}|$ is small (i.e. the viscosity gradient is small, as at late times), ${\mathcal{S}}$ becomes independent of $\unicode[STIX]{x1D6FD}$ and tends to the value for $R=0$ . This suggests that at late enough times, the effective diffusivity will always become independent of the log-viscosity ratio, and the flow will always evolve like the neutrally stable case.

Figure 7. Plot of ${\mathcal{S}}$ versus $|\unicode[STIX]{x1D6FD}|$ for $R<0$ and $R>0$ . The solid black lines correspond to ${\mathcal{S}}$ , equation (4.8), calculated by numerically integrating (4.7). The leading-order small- $|\unicode[STIX]{x1D6FD}|$ asymptotic behaviour, which is equal to when $R=0$ , is given by the dashed red line and the next-order corrections in $\unicode[STIX]{x1D6FD}$ are given by the dot-dashed blue lines, equation (4.12).

When $|R|<\unicode[STIX]{x1D70E}$ , we need to only consider the small- $|\unicode[STIX]{x1D6FD}|$ limit because by the time the flow reaches the late-time regime, $\unicode[STIX]{x1D6FD}$ is inevitably small. We can justify this claim using a scaling argument: once the flow reaches the late-time regime, which occurs at a time $t=O(Pe)$ , the mixing length will have grown to a width $h=O(Pe\unicode[STIX]{x0394}u)$ , or equivalently, the concentration gradient is $\unicode[STIX]{x2202}\overline{c}/\unicode[STIX]{x2202}x=O(1/Pe\unicode[STIX]{x0394}u)$ . This means that when the flow transitions to the late-time regime, $|\unicode[STIX]{x1D6FD}|=O(R/\unicode[STIX]{x0394}u)<O(1)$ .

4.3.1 Small viscosity gradient limit: $|\unicode[STIX]{x1D6FD}|\ll 1$

For $\unicode[STIX]{x1D6FD}\ll 1$ , we start by expanding the concentration deviations as $\tilde{c}=\unicode[STIX]{x1D6FD}\tilde{c}_{1}+\unicode[STIX]{x1D6FD}^{2}\tilde{c}_{2}+O(\unicode[STIX]{x1D6FD}^{3})$ . Substituting into (4.6), expanding and equating powers of $\unicode[STIX]{x1D6FD}$ gives

(4.9a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\tilde{c}_{1}}{\unicode[STIX]{x2202}y^{2}}=k-1,\quad \frac{\unicode[STIX]{x2202}^{2}\tilde{c}_{2}}{\unicode[STIX]{x2202}y^{2}}=(k\tilde{c}_{1}-k\overline{k\tilde{c}_{1}}).\end{eqnarray}$$

Given that the concentration deviations must satisfy periodicity and have vanishing mean, we find that

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{c}_{1}={\mathcal{I}}(k-1)-\overline{{\mathcal{I}}}(k-1), & \displaystyle\end{eqnarray}$$

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{c}_{2}={\mathcal{I}}(k\tilde{c}_{1})-\overline{{\mathcal{I}}}(k\tilde{c}_{1})-\overline{k\tilde{c}_{1}}{\mathcal{I}}(k)+\overline{k\tilde{c}_{1}}\,\overline{{\mathcal{I}}}(k), & \displaystyle\end{eqnarray}$$

where ${\mathcal{I}}(f)=\int _{0}^{y}\int _{0}^{\unicode[STIX]{x1D701}}f\,\text{d}\unicode[STIX]{x1D702}\,\text{d}\unicode[STIX]{x1D701}$ for a given function $f$ and, as before, the overbar refers to a transverse average. The shear-enhanced dispersivity, ${\mathcal{S}}$ , is thus

(4.12) $$\begin{eqnarray}{\mathcal{S}}=-\frac{1}{\unicode[STIX]{x1D6FD}}\int _{0}^{1}u\tilde{c}\,\text{d}y=\overline{k\tilde{c}_{1}}+\unicode[STIX]{x1D6FD}(\overline{k\tilde{c}_{1}^{2}}-\overline{k\tilde{c}_{1}}^{2}+\overline{k\tilde{c}_{2}})+O(\unicode[STIX]{x1D6FD}^{2}).\end{eqnarray}$$

The leading-order contribution to ${\mathcal{S}}$ is identical to the $R=0$ limit in § 3.3 and is plotted in figure 7 as a dashed red line. The first-order corrections are given by dot-dashed blue lines.

4.3.2 Comparison to numerical simulations

We now use the calculated shear-enhanced dispersivity to determine how $\overline{c}$ evolves in time. Allowing ${\mathcal{S}}$ to vary in space, (2.16) yields a nonlinear diffusion equation for $\overline{c}$ ,

(4.13) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}t}=\frac{1}{Pe}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left\{[1+Pe^{2}{\mathcal{S}}(\unicode[STIX]{x1D6FD}(x))]\frac{\unicode[STIX]{x2202}\overline{c}}{\unicode[STIX]{x2202}x}\right\},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}(x)=R\,Pe\,\unicode[STIX]{x2202}\overline{c}/\unicode[STIX]{x2202}x$ depends on the local mean concentration gradient. We solve (4.13) numerically with no-flux boundary conditions in the far field using a Crank–Nicolson predictor–corrector method. We use both the exact form for ${\mathcal{S}}$ (by solving (4.7) and calculating (4.8)), and the small- $|\unicode[STIX]{x1D6FD}|$ approximation in (4.12). The model concentration fields are initialized with a diffuse error-function solution although the long-time results are indifferent to the exact initial conditions. The non-uniform dispersion results in profiles that deviate from the classical error-function solution owing to the enhanced or diminished dispersion in regions of large concentration gradients.

Figure 8. Evolution of the mixing length and transversely averaged concentration for $|R|<\unicode[STIX]{x1D70E}$ in the late-time regime. (a) Evolution of the mixing length $h$ , normalized by the neutral displacement mixing length $h_{R=0}$ for $(\unicode[STIX]{x1D70E},Pe,n)=(1,100,1)$ and $R$ ranging from $-0.4$ to $0.4$ . (b) Plot of the difference between $\overline{c}$ for $R\neq 0$ and $R=0$ . The profiles are measured at $t=100$ for the same parameters as (a). The theoretical predictions, found by solving (4.13) with either the exact solution of ${\mathcal{S}}$ (by solving (4.7)), or (4.12), are given by the dashed and dotted black lines respectively.

Figure 8(a) shows the evolution of the normalized mixing length, $h/h_{R_{0}}$ , where $h_{R_{0}}=h(t,R=0)$ , for different $R$ , from the full two-dimensional (2-D) numerical simulations (solid coloured lines) and model results (dashed and dotted black lines). As expected, increasing $R$ results in increased spreading, while the effects of variations in the viscosity reduce as $t$ is increased. The reduced-order model not only captures this behaviour but also accurately predicts the manner in which the flow evolves. This can be further seen in figure 8(b) which compares the reduced-order model predictions for $\overline{c}$ with the full 2-D numerical simulations. The very good agreement between both of the model solutions and the numerical results suggests that the late-time behaviour can be accurately modelled by (4.13) with (4.12).

5.1 Concentration model

As found in § 4.3, the shear-enhanced dispersivity ${\mathcal{S}}$ is given by (4.8) where $\unicode[STIX]{x1D6FD}=R\,Pe\,\unicode[STIX]{x2202}\overline{c}/\unicode[STIX]{x2202}x$ and $\tilde{c}$ is given by (4.7) and $u$ by (4.3). However, unlike in § 4.3, for viscosity-stabilized flows, we now expect this model to apply for all time, rather than just at late times, given that the flow remains nearly transversely uniform ( $c^{\prime }\ll 1$ ). In particular, $|\unicode[STIX]{x1D6FD}|$ is not just small but can take on any value. We thus first consider the large- $\unicode[STIX]{x1D6FD}$ limit.

We start by expanding $\tilde{c}$ in powers of $1/\unicode[STIX]{x1D6FD}$ , $\tilde{c}=\tilde{c}_{0}+\tilde{c}_{1}/\unicode[STIX]{x1D6FD}+O(1/\unicode[STIX]{x1D6FD}^{2})$ . Substituting into (4.7) expanding and equating different powers of $\unicode[STIX]{x1D6FD}$ , we find

(5.1) $$\begin{eqnarray}\frac{k\text{e}^{\tilde{c}_{0}}}{\displaystyle \int _{0}^{1}k\text{e}^{\tilde{c}_{0}}\,\text{d}y}-1=0,\end{eqnarray}$$

and so,

(5.2a,b ) $$\begin{eqnarray}\tilde{c}_{0}=-\ln (k)+\overline{\ln (k)},\quad \tilde{c}_{1}=\frac{\unicode[STIX]{x2202}^{2}\tilde{c}_{0}}{\unicode[STIX]{x2202}y^{2}},\end{eqnarray}$$

and

(5.3) $$\begin{eqnarray}\tilde{c}_{1}=-\left(\frac{kk^{\prime \prime }-(k^{\prime })^{2}}{k^{2}}\right),\end{eqnarray}$$

where we have again used the fact that $\int _{0}^{1}\tilde{c}_{0}=\int _{0}^{1}\tilde{c}_{1}=0$ . To leading order, the concentration deviations align themselves such that the streamwise velocity is zero. The shear-enhanced dispersivity in this limit is

(5.4) $$\begin{eqnarray}{\mathcal{S}}=\frac{1}{\unicode[STIX]{x1D6FD}}\int _{0}^{1}u\tilde{c}\,\text{d}y=\frac{1}{\unicode[STIX]{x1D6FD}^{2}}\int _{0}^{1}\ln (k)\left(\frac{kk^{\prime \prime }-(k^{\prime })^{2}}{k^{2}}\right)\,\text{d}y+O\left(\frac{1}{\unicode[STIX]{x1D6FD}^{3}}\right).\end{eqnarray}$$

In the case of sinusoidally varying log-permeability, this limit corresponds to

(5.5) $$\begin{eqnarray}{\mathcal{S}}=\frac{2\unicode[STIX]{x1D70E}^{2}\unicode[STIX]{x03C0}^{2}}{\unicode[STIX]{x1D6FD}^{2}}+O\left(\frac{1}{\unicode[STIX]{x1D6FD}^{3}}\right).\end{eqnarray}$$

Figure 10. (a) Plot of ${\mathcal{S}}$ normalized by its small- $|\unicode[STIX]{x1D6FD}|$ limit versus $\unicode[STIX]{x1D6FD}$ and (b) plot of ${\mathcal{S}}$ normalized by its large- $|\unicode[STIX]{x1D6FD}|$ limit versus $\unicode[STIX]{x1D6FD}$ for $\unicode[STIX]{x1D70E}$ ranging from $1$ to $5$ . The dashed black lines denote the asymptotic limits (4.12) and (5.4) in (a) and (b) respectively and the dotted lines correspond to the approximate solution (5.6) for $\unicode[STIX]{x1D70E}=5$ .

Equations (4.12) and (5.4) correspond to the asymptotic limits of small and large viscosity gradient, and are plotted together with the full solutions of (4.7)–(4.8) for the shear-enhanced dispersivity in figures 10(a) and 10(b) respectively. We can also combine the two limits into a very simple approximate analytical composite solution; for example for the case of a sinusoidally varying log-permeability,

(5.6) $$\begin{eqnarray}{\mathcal{S}}_{comp}=\left(\frac{1}{{\mathcal{S}}^{\ast }}+\frac{\unicode[STIX]{x1D6FD}^{2}}{2\unicode[STIX]{x1D70E}^{2}\unicode[STIX]{x03C0}^{2}}\right)^{-1},\end{eqnarray}$$

where ${\mathcal{S}}^{\ast }=\int _{0}^{1}[\int _{0}^{y}(k-1)\,\text{d}\unicode[STIX]{x1D702}]^{2}\,\text{d}y$ is the leading-order behaviour in (4.12). This approximate solution captures the general behaviour of the shear-enhanced dispersivity reasonably well without having to solve (4.7) exactly (see dotted lines in figure 10).

5.2 Comparison to numerical simulations

We solve the reduced-order model (4.13) both directly and using the approximate composite solution (5.6), in the same manner as § 4.3.2 with a step initial condition (2.15). The comparisons to the transversely averaged concentration profiles from the full numerical simulations are given in figure 11(b). The exact solution of the reduced model is not only able to capture the sharp interface and the long tails, but it also accurately predicts the evolution of the concentration field. Using the composite approximation of ${\mathcal{S}}$ in the reduced model also captures the qualitative evolution of the transversely averaged concentration although it slightly overestimates the amount of spreading that occurs. The mixing length is also shown as a function of time for different $\unicode[STIX]{x1D70E}$ in figure 11(a). The very good agreement between the reduced-order model and the full simulations over a range of $\unicode[STIX]{x1D70E}$ and $t$ , suggests that the full 2-D problem can be reduced to solving a 1-D nonlinear diffusion equation for all times in this limit.

Figure 11. Evolution of the mixing length and transversely averaged concentration for $|R|>\unicode[STIX]{x1D70E}$ , $R<0$ . (a) Plot of $h$ versus $t$ for $(R,Pe,n)=(-2,300,1)$ and $\unicode[STIX]{x1D70E}$ ranging from $0$ to $0.3$ . (b) Plot of $\overline{c}(x,t)$ for $(R,\unicode[STIX]{x1D70E},Pe,n)=(-2,0.2,300,1)$ and $t=32,100,320,1000$ . The theoretical predictions, found by solving (4.13), with either the exact solution of ${\mathcal{S}}$ (found by solving (4.7)), or (5.6), are given by the dashed and dot-dashed black lines respectively.

6.1 Regime I: fingering within layers

The flow is able to finger within the permeability layers when the length scale of the most unstable mode, $y\sim 1/RPe$ (Tan & Homsy Reference Tan and Homsy1986), is small compared to the length scale imposed by the permeability $y\sim 1/n$ . Hence the flow always fingers for sufficiently large $Pe$ and is also observed to be enhanced when $\unicode[STIX]{x1D70E}$ is small (De Wit & Homsy Reference De Wit and Homsy1997b ; Shahnazari et al. Reference Shahnazari, Maleka Ashtiani and Saberi2018). Note that fingering within layers is also possible when permeability effects dominate over viscous effects, $R<\unicode[STIX]{x1D70E}$ , so long as $Pe$ is sufficiently large.

Figure 13. Evolution of the concentration field for $|R|>\unicode[STIX]{x1D70E}$ and $R>0$ in the fingering within layers regime; $(R,Pe,n)=(2,8000,2)$ . (a,b) Colour maps of the concentration field for (a) $\unicode[STIX]{x1D70E}=0.1$ and (b) $\unicode[STIX]{x1D70E}=0$ at $t=1$ . (c) Plot of the transversely averaged concentration $\overline{c}(x,t)$ and (d) transverse variance in concentration $var(c)=\int _{0}^{1}{c^{\prime }}^{2}\,\text{d}y$ versus the similarity variable $x/t$ for $\unicode[STIX]{x1D70E}$ ranging from $0$ to $0.3$ . In (c), the theoretical prediction for a homogeneous medium (Nijjer et al. Reference Nijjer, Hewitt and Neufeld2018, dashed) and for a heterogeneous medium ((C 4), dot-dashed) are also shown.

Figure 13(a,b) shows snapshots of the concentration field for $\unicode[STIX]{x1D70E}=0.1$ and $\unicode[STIX]{x1D70E}=0$ respectively. The fingers evolve in a similar manner in both cases, leading to an asymmetric transversely averaged concentration that evolves self-similarly with similarity variable $x/t$ (figure 13 c). The main difference between the two simulations is that the fingers move significantly faster when permeability heterogeneities are present. For example, the mixing length grows at more than double the rate in the simulation in figure 13(a) compared to the homogeneous case, which is far larger than one would predict simply by considering the difference in permeability.

This difference is due to the structure of the flow being fundamentally different in the two simulations. Whereas, in the homogeneous medium, the forward and backward propagating fingers are on average about the same size, in the heterogeneous medium, the backward propagating fingers are broad and aligned with the low permeability layers. This means the two fluids tend to be much more segregated (figure 13 d), and so the effective viscosity difference between the forward- and backward-propagating fingers is much larger in the heterogeneous medium, which leads to much faster spreading. We also observe from the simulations that the spreading rate is broadly insensitive to $\unicode[STIX]{x1D70E}$ for $\unicode[STIX]{x1D70E}\gtrsim 0.1$ (figure 13 c,d), which supports the idea that the presence of layered heterogeneity changes the spreading rate through the structure of the flow. In § C.1 we extend the model presented by Nijjer et al. (Reference Nijjer, Hewitt and Neufeld2018) for the homogeneous case by accounting for this change in finger structure. The model prediction, given by the dot-dashed line in figure 13(c), gives good quantitative agreement with the numerical simulations when $\unicode[STIX]{x1D70E}\gtrsim 0.1$ . Thus we find that, although small amounts of permeability heterogeneity would be expected to have a small effect on the spreading rate, they can in fact alter the structure of the flow and lead to significantly faster spreading and mixing.

6.2 Regime II: stable channelling

After the fluid fingers inside the high-permeability layers, it coarsens until a single broad finger remains in each layer (figure 12 b). Figure 14(a–d) shows a sequence of snapshots of the concentration field with overlain streamlines for $(R,\unicode[STIX]{x1D70E},Pe,n)=(2,0.1,1000,1)$ for which at intermediate times the flow is always in this channelling regime. In this example, $n=1$ and so the channelling consists of a single pair counter-propagating fingers. The flow is predominantly longitudinal except in regions localized near large concentration gradients. Initially the streamwise velocity is large and localized in the finger. Upstream and downstream, where the viscosity is uniform, the velocity is imposed by the permeability, but is small relative to the contributions due to viscosity variations (figure 14 f). As time progresses, and the fluids become more mixed, the effect of the viscosity reduces and the streamwise velocity in the fingered region approaches the upstream and downstream velocity (figure 14 b–d). The concentration deviations, figure 14(e), are uniform in the $x$ -direction, vary sinusoidally in the $y$ direction and are in phase with the velocity.

Figure 14. Evolution of the concentration field for $|R|>\unicode[STIX]{x1D70E}$ and $R>0$ in the channelling regime; $(R,\unicode[STIX]{x1D70E},Pe,n)=(2,0.1,500,1)$ . (a–d) Colour maps of the concentration field with overlain streamlines at (a) $t=20$ , (b) $t=60$ , (c) $t=100$ and (d) $t=140$ . (e) Colour map of the concentration deviations $c^{\prime }=c-\overline{c}$ at $t=20$ . (f) Colour map of the streamwise velocity $u$ at $t=20$ . Note that the aspect ratio of the figures is compressed by a factor of 30, so variations in the $x$ -direction seem more pronounced than they actually are.

Figure 15. (a) Plot of $\overline{c}(x,t)$ from the numerical simulations (solid lines) and the theoretical prediction overlain (dashed line) for $(R,\unicode[STIX]{x1D70E},Pe,n)=(2,0.1,1000,1)$ . (b) Plot of the linear slope $\unicode[STIX]{x1D6FC}$ of $\overline{c}$ , as a function of $\unicode[STIX]{x1D70E}$ , showing a least-squares fit to the data from numerical simulations (dots) and the theoretical solution (C 11) for $K=0.5$ (solid line) and $K=0.6$ (dashed line).

This behaviour resembles the late-time regime in the miscible viscous-fingering instability (i.e. when $\unicode[STIX]{x1D70E}=0$ ). More specifically, the dynamics consists of an interior region where the background gradient is linear and steady and the ends of which are filled in by the propagating fingers (figure 15 a). Superimposed on the steady background gradient are sinusoidal concentration deviations which decay in time.

In § C.2, we generalize the analysis from § 4.3.1 by allowing the deviations to evolve in time, and so derive a simplified model for the evolution of the mean concentration field in this regime. To test the validity of this model, we compare the predicted steady interior slope, $\unicode[STIX]{x1D6FC}$ , to the results of the 2-D numerical simulations (figure 15 b). As the viscosity ratio and permeability variance are increased $\unicode[STIX]{x1D6FC}$ decreases due to the increased fluid velocity. The model shows reasonable agreement with the numerical simulations and the slight over-prediction of the slope is likely due to the underestimation of the velocity at early times.

6.3 Regimes III and IV: viscous fingering across layers

As noted earlier, if $n>1$ , the channelling regime can become unstable to a viscous-fingering instability which has fingers that are wider than the permeability structure imposed. The viscous fingers that develop in the heterogeneous medium are qualitatively similar to the fingers that develop in a homogeneous medium (figure 16 a,b) and go through the same large-scale coarsening until a single finger remains. This single finger then evolves in manner similar to the single-finger state in a homogeneous medium ( $\unicode[STIX]{x1D70E}=0$ ).

Figure 16. Evolution of the concentration field for $|R|>\unicode[STIX]{x1D70E}$ and $R>0$ in the viscous fingering across layers regime; $(R,Pe,n)=(2,2000,32)$ . (a,b) Colour maps of the concentration field for (a) $\unicode[STIX]{x1D70E}=0.1$ and (b) $\unicode[STIX]{x1D70E}=0$ at $t=8$ . (c) Plot of $h(t)$ for five different simulations each (light, thin lines) and their average (dark, thick line). (d) Plot of $\overline{c}(x,t)$ versus the similarity variable $x/t$ and $t$ ranging from 6 to 20. Each curve corresponds to the average across five different simulations.

One key difference between the two simulations is that in the heterogeneous medium there tends to be more tip splitting, aligned with the permeability field, due to velocity heterogeneities at the finger tips as well as fading and coalescence of fingers. The fact that the fingers are more unstable leads to more variability from simulation to simulation, more intermittent flow and non-uniform growth of the mixing length (figure 16 c). The concentration field evolves in qualitatively the same manner in both cases, the profile is asymmetric and again has the similarity variable $x/t$ (figure 16 d). However the spreading rate is slower in the case when heterogeneities are present. Modelling the difference in spreading rate is the subject of future work.

Note that viscous fingering across layers is also possible when permeability effects dominate over viscous effects, $R<\unicode[STIX]{x1D70E}$ . This can occur when the system reaches the late-time regime and the interface is still hydrodynamically unstable. For the interface to be stable, the width of the mixing zone must be $h\sim RPe$ (Nijjer et al. Reference Nijjer, Hewitt and Neufeld2018). The mixing zone at the start of the late-time regime in the permeability-dominated scenario is $h\sim Pe\unicode[STIX]{x0394}u/n^{2}\sim Pe\unicode[STIX]{x1D70E}/n^{2}$ and so the late-time state may be unstable to further coarsening via a viscous instability if $R>\unicode[STIX]{x1D70E}^{2}/n^{2}$ .