Stable and unstable miscible displacements in layered porous media

Abstract

The effect of permeability heterogeneities and viscosity variations on miscible displacement processes in porous media is examined using high-resolution numerical simulations and reduced theoretical modelling. The planar injection of one fluid into a fluid-saturated, two-dimensional porous medium with a permeability that varies perpendicular to the flow direction is studied. Three cases are considered, in which the injected fluid is equally viscous, more viscous or less viscous than the ambient fluid. In general it is found that the flow in each case evolves through three regimes. At early times, the flow exhibits the concentration evolves diffusively, independent of both the permeability structure and the viscosity ratio. At intermediate times, the flow exhibits different dynamics including channelling and fingering, depending on whether the injected fluid is more or less viscous than the ambient fluid, and depending on the relative magnitude of the viscosity and permeability variations. Finally, at late times, the flow becomes independent of the viscosity ratio and dominated by shear-enhanced (Taylor) dispersion. For each of the regimes identified above, we develop reduced-order models for the evolution of the transversely averaged concentration and compare them to the full numerical simulations.

1 Introduction

Miscible displacement processes, in which a viscous fluid is injected into a fluid-saturated porous medium, are relevant in a number of problems including enhanced oil recovery (Lake 1989), carbon capture and storage (Huppert & Neufeld 2014) and subsurface contaminant transport (Abriola 1987). In many of these cases, the goal is to optimize and control the displacement front and the amount of mixing that occurs between the two fluids. The evolution of the displacement front is controlled by a number of factors; in this paper we will examine the combined effects of both permeability heterogeneities and viscosity variations.

If the porous medium is homogeneous and the injected and ambient fluids have different viscosities then, depending on the viscosity ratio, the displacement front can be stable or unstable. If the injected fluid is more viscous, the interface is stable and the fluid front moves uniformly, whereas if the injected fluid is less viscous the displacement front can be unstable and lead to complex fingering patterns. Miscible viscous fingering has been extensively studied using a number of different approaches (for example Tan & Homsy 1986, 1987, 1988; Zimmerman & Homsy 1991, 1992; Jha, Cueto-Felgueroso & Juanes 2011a , b ; Chui, De Anna & Juanes 2015; Pramanik & Mishra 2015b ). Recently, Nijjer, Hewitt & Neufeld (2018) described the full lifecycle of miscible viscous fingering: they identified the different regimes through which the flow evolves and in each regime modelled the evolution of the flow. We take a similar approach in this paper and look at the role of permeability heterogeneities on the temporal evolution of the displacement front in the case of neutral, stable and unstable displacements.

Most physically relevant porous media are not homogeneous but vary on a wide range of length scales from the pore scale to the reservoir scale in both an ordered and disordered manner (Weber 1986). A number of studies have considered the combined effects of both randomly varying permeability fields and viscosity variations on miscible displacements using theoretical (Welty & Gelhar 1991), numerical (Tan & Homsy 1992; Waggoner, Castillo & Lake 1992; Chen & Meiburg 1998; Camhi, Meiburg & Ruith 2000; Talon et al. 2004; Tchelepi et al. 2004; Nicolaides et al. 2015) and experimental (Jiao & Hotzl 2004; Tchelepi et al. 2004) approaches. These studies have demonstrated that the flow exhibits a range of dynamical behaviour, including dispersing, channelling and fingering, and highlight the complexity of the flow patterns that arise.

In this work, we focus on large-scale permeability variations that are perpendicular to the flow direction. This structure is widespread in nature, being characteristic of geological formations consisting of different sedimentary sequences. Previous work looking at the combined effect of permeability variations and an unstable flow configuration found that adding permeability heterogeneities tends to cause channelling; that is, flow predominantly along permeability layers rather than chaotic fingering (De Wit & Homsy 1997b ; Sajjadi & Azaiez 2013; Shahnazari, Maleka Ashtiani & Saberi 2018). Loggia et al. (1996) found that when the injected fluid is more viscous than the ambient channelling is observed when the viscosity ratio is smaller than the ratio of permeabilities, and a single dispersive front is attained when the viscosity ratio is larger than the ratio of permeabilities. While these studies highlight some of the interesting qualitative behaviour that can be observed in miscible displacement flows when heterogeneity and viscosity variations interact, they do not provide a full overview of the different dynamical regimes that occur and the temporal evolution of the flow between them. The aim of this work is to identify the full lifecycle and evolution of stable and unstable miscible displacements in layered porous media, and to develop reduced-order models for the spreading and dispersion of the fluids, which can be used to quantitatively predict and up-scale flow in heterogeneous porous media.

This paper is laid out as follows. In § 2, we formulate the problem and outline the numerical method we use to solve it. In § 3, we consider neutrally stable displacements where the two fluids have the same viscosity. In § 4 we consider the effect of small viscosity variations, both stabilizing and de-stabilizing. In § 5 we consider large stabilizing viscosity variations and in § 6 we consider large de-stabilizing viscosity variations. In each of §§ 3–6 we discuss the time evolution of the concentration field and the different regimes through which the flow evolves, and derive reduced-order models for the evolution of the concentration field.

2 Problem formulation

A schematic of the problem geometry is given in figure 1. We consider a semi-infinite, two-dimensional porous strip with finite width $a$ . For simplicity, the porosity of the medium, $\unicode[STIX]{x1D719}$ , is constant, but the permeability, $k=k(y)$ , varies in the direction perpendicular to the flow. A fluid of viscosity $\unicode[STIX]{x1D707}_{1}$ is injected into the medium at a constant volumetric flux $Q$ , which is initially saturated with an ambient fluid of viscosity $\unicode[STIX]{x1D707}_{2}$ . The two fluids are fully miscible.

Figure 1. A schematic of the model geometry. The porous medium is semi-infinite with width $a$ and has a permeability structure that is only a function of the transverse coordinate, $y$ . The porous medium is initially saturated with a fluid of viscosity $\unicode[STIX]{x1D707}_{2}$ . Another fluid, with viscosity $\unicode[STIX]{x1D707}_{1}$ , which is fully miscible with the first, is injected at a constant flux $Q$ .

We assume that the flow obeys Darcy’s law everywhere, such that the pore-scale Reynolds number is small, the flow is incompressible, and the concentration, $c$ , of the injected fluid evolves via advection and diffusion;

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=-\frac{k(y)}{\unicode[STIX]{x1D707}(c)}\unicode[STIX]{x1D735}p, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c=\unicode[STIX]{x1D719}D\unicode[STIX]{x1D6FB}^{2}c. & \displaystyle\end{eqnarray}$$

Here $\boldsymbol{u}=(u,v)$ is the Darcy velocity and $p$ is the pressure, and $D$ is the effective diffusivity or dispersion coefficient between the two fluids, which in porous media tends to be velocity dependent, but is here assumed, for simplicity, to be a constant. The viscosity varies with the concentration of the injected fluid, and, following previous work (e.g. Tan & Homsy 1986), we assume an Arrhenius-like dependence,

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D707}(c)=\unicode[STIX]{x1D707}_{2}\text{e}^{-Rc},\end{eqnarray}$$

where $R=\log (\unicode[STIX]{x1D707}_{2}/\unicode[STIX]{x1D707}_{1})$ .

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 1997a , b ; Sajjadi & Azaiez 2013), 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 1986; Nijjer et al. 2018), 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. 2011a ), 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. (2018). 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 (1999). 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 (2015a ), 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 Neutrally stable displacements $R=0$

We first consider the case where the injected and ambient fluids have the same viscosity. While this case has been explored by a number of authors (e.g. Camacho 1993; Berentsen, Verlaan & van Kruijsdijk 2005; Dentz & Carrera 2007), we outline it here both for completeness and to set the stage for the analysis in the remainder of the manuscript. Results from a representative simulation are given in figure 2. As was noted earlier, if the flow is hydrodynamically stable, $n$ can be scaled out of the problem and so we only consider the case where $n=1$ . In this case, the permeability is highest in the centre of the channel and lowest at the top and bottom boundaries (figure 2 a). In the moving frame, this causes the fluid in the middle of the channel to move to the right and the fluid near the top and bottom boundaries to move to the left. This shear spreads and mixes the fluids. To quantify this spreading, we plot the evolution of the mixing length, $h(t)$ , in figure 2(b). We find that the mixing evolves through three distinct regimes each with a different scaling behaviour. The concentration fields corresponding to each of these regimes are plotted in figure 2(c–e). In the first and third regimes, the concentration fields look nearly indistinguishable: the concentration is nearly transversely uniform and relatively diffuse in the streamwise direction. In the second regime, there is a relatively sharp interface aligned with the permeability variations. Based on these observations, we expect that in the first (early-time) and third (late-time) regimes, the streamwise transport is diffusively dominated, whereas in the second (intermediate-time) regime, the streamwise transport is advectively dominated.

Since $R=0$ , the concentration acts as a passive tracer. This means that the velocity is decoupled from the concentration field, and is given from (2.5) by

(3.1a,b ) $$\begin{eqnarray}u(y)=k(y)-1,\quad v=0.\end{eqnarray}$$

In the case of sinusoidally varying log-permeability (2.9), the velocity is,

(3.2a,b ) $$\begin{eqnarray}u=\frac{\text{e}^{-\unicode[STIX]{x1D70E}\cos (2\unicode[STIX]{x03C0}y)}}{\text{I}_{0}(\unicode[STIX]{x1D70E})}-1,\quad v=0.\end{eqnarray}$$

Given this fixed, known velocity the concentration simply evolves via the advection–diffusion equation (2.7). In the following sections, we consider the dominant balances in (2.7) to determine how the concentration field evolves in time.

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 1953; Aris 1956). 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 1983).

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 Small viscosity variations $|R|<\unicode[STIX]{x1D70E}$

Next we consider the effect of viscosity variations that are weak compared to the permeability; that is, the log-viscosity ratio is smaller than the log-permeability variance, $|R|<\unicode[STIX]{x1D70E}$ .

Figure 5. Colour maps of the concentration field with overlain streamlines for $|R|<\unicode[STIX]{x1D70E}$ . (a,b) $R=0.4$ , (c,d) $R=0$ and (e,f) $R=-0.4$ and $(\unicode[STIX]{x1D70E},Pe,n)=(1,500,1)$ . The snapshots are taken at (a,c,e) intermediate times ( $t=1$ ) and (b,d,f) late times ( $t=31$ ). Note that the aspect ratios of the late-time figures are distorted.

In the absence of permeability variations, when $R>0$ and the Péclet number is sufficiently large, the flow is unstable and a set of complex nonlinearly evolving fingers develop (Tan & Homsy 1988). If permeability layering is introduced, the flow tends to be forced along the permeability pathways (De Wit & Homsy 1997b ; Shahnazari et al. 2018) and as the permeability variance is increased, the flow becomes more and more channelized until the fingers no longer interact. This is especially true when the permeability variability dominates over variations in viscosity, $\unicode[STIX]{x1D70E}\gg |R|$ . Although instabilities are still possible (and are further discussed in § 6), in this section we focus on flows that remain hydrodynamically stable and follow the permeability pathways imposed.

Figure 5 shows the concentration field overlain with streamlines for $\unicode[STIX]{x1D70E}=1$ and $R=0.4,0$ and $-0.4$ at intermediate times (left) and late times (right). For all three values of $R$ , we find that the concentration field evolves in qualitatively the same manner: after an early-time diffusive regime, as in § 3.1, at intermediate times the flow is dominated by advective stretching (figure 5 a,c,e); and at late times the flow is dominated by shear-enhanced dispersion (figure 5 b,d,f). The main difference between flows where $R\neq 0$ and $R=0$ is that at intermediate times the interface is either stretched ( $R>0$ ) or compressed ( $R<0$ ) relative to the neutrally stable case owing to the viscosity-enhanced or viscosity-tempered streamwise velocity. At late times, the viscosity contrast seems to have little effect and the concentration field and streamlines look nearly indistinguishable. In the following subsections we examine the effects of small viscosity variations on the evolution of the three regimes identified in § 3. We begin in § 4.1 with a consideration of how viscosity contrasts affect the fluid velocity.

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 1995), 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 Large stabilizing viscosity variations $|R|>\unicode[STIX]{x1D70E},R<0$

Whereas in the previous section the flow evolves in qualitatively the same manner as the uniform viscosity case ( $R=0$ ), when the viscosity ratio is larger than the permeability ratio, the concentration evolves in a qualitatively different manner. Here we consider the case where the injected fluid is more viscous than the ambient fluid ( $R<0$ ) and the magnitude of the log-viscosity ratio is larger than the log-permeability ratio.

Figure 9. Evolution of the concentration field for stable displacements $|R|>\unicode[STIX]{x1D70E}$ , $R<0$ . $(R,\unicode[STIX]{x1D70E},Pe,n)=(-1,0.1,100,1)$ . (a–d) Colour maps of the concentration field with overlain streamlines at (a) $t=3.6$ , (b) $t=17.8$ , (c) $t=89.1$ and (d) $t=447$ . (e) Colour map of the concentration deviations $c^{\prime }=c-\overline{c}$ at $t=3.6$ . (f) Colour map of the streamwise velocity $u$ at $t=3.6$ .

Figure 9(a–d) shows a sequence of snapshots of the concentration field in this limit. The interface is initially not stretched by the permeability variations owing to the large longitudinal gradient in viscosity (figure 9 a). This is because, for large $|\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}x|$ , distorting the interface generates large transverse gradients in concentration. These correspond to large, stable transverse gradients in viscosity (and hence pressure) which tend to force the profile back to vertical. Ultimately, this results in a stationary interface (in the moving frame) where the streamwise velocity goes to zero (figure 9 f). Far upstream and downstream, where the viscosity is transversely uniform, the velocity is imposed by the permeability. The abrupt change in velocity at the interface drives a circulation on either side of the interface carrying concentration away from it. As the fluids mix, the concentration, and hence viscosity gradients at the interface weaken. Over time, the stabilizing viscous forces become sufficiently weak that the streamwise velocity can grow and streamlines begin to penetrate through the interface (figure 9 b,c). Eventually the interface becomes very diffuse and the velocity becomes predominantly longitudinal (figure 9 d) as in the cases considered in previous sections. In contrast to those cases, however, where the concentration deviations, ( $c^{\prime }$ ), are $O(1)$ until the late-time regime, the concentration deviations here are always small because the concentration remains relatively transversely uniform with no large-scale channelling into layers.

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 Large de-stabilizing viscosity variations $|R|>\unicode[STIX]{x1D70E},R>0$

Finally, we consider the case where the injected fluid is less viscous than the ambient fluid and the log-viscosity ratio is larger than the log-permeability ratio, $R>\unicode[STIX]{x1D70E}$ . In the absence of any permeability variations, $\unicode[STIX]{x1D70E}=0$ , this configuration is hydrodynamically unstable (provided $Pe$ is sufficiently large). In homogeneous media such unstable miscible displacements evolve through three flow regimes (Nijjer et al. 2018): at early times, the flow is linearly unstable, where the interface grows diffusively, and fingers grow exponentially; at intermediate times, the fingers have finite amplitude and propagate and interact with each other nonlinearly leading to coarsening; and at late times, a single pair of counter-propagating fingers remain which propagate and slow, leaving a well-mixed interior.

Figure 12. Evolution of the concentration field for $|R|>\unicode[STIX]{x1D70E}$ and $R>0$ . (a–d) Colour maps of the concentration field for $(R,\unicode[STIX]{x1D70E},Pe,n)=(2,0.1,4000,4)$ at (a) $t=0.1$ , (b) $t=8$ , (c) $t=104$ and (d) $t=1500$ . Note that the aspect ratio is increasingly compressed in figures (b–d). (e,f) Evolution of $h(t)$ for $(Pe,n)=(1000,4)$ and: (e) $R=2$ with varying $\unicode[STIX]{x1D70E}$ ; and (f) $\unicode[STIX]{x1D70E}=0.1$ with varying $R$ .

In the presence of permeability layering, there is a competition between the evolving wavelength of viscous fingering and the imposed wavelength of the permeability structure. The competition between viscous fingering, which acts to coarsen the transverse length scale, and permeability layering, which acts to impose a fixed length scale, results in rich intermediate-time dynamics. As a result, the number of layers $n$ can no longer be scaled out of the problem. Nonetheless, as before, the early-time dynamics is still dominated by longitudinal diffusion across the sharp interface, and the late-time dynamics is dominated by shear-enhanced dispersion which becomes independent of the viscosity ratio at long times (cf. § 4.3).

In general, there are four possible intermediate-time regimes through which the flow can evolve, representative snapshots of which are given in figure 12. In the first regime (I), fingering occurs within the permeability layers and these fingers coarsen until they coincide with the imposed permeability layering (figure 12 a). In the second regime (II), the flow follows the imposed layered structure while diffusing across the layers causing the flow to slow down (figure 12 b). In some cases, this flow can then be unstable, leading to a third regime (III) which corresponds to fingering over a transverse length scale that is larger than that imposed by the permeability (figure 12 c). These fingers then also coarsen, leading to a fourth regime (IV) where a single pair of counter-propagating fingers remain (figure 12 d). This pair of fingers slows leaving a well-mixed region which eventually evolves through shear-enhanced dispersion. Note that regimes III and IV can only occur if $n>1$ .

The growth of the mixing length can be drastically different depending on which regime the flow is in (figure 12 e,f). When the flow fingers (regimes I and III), as with a homogeneous porous medium, the mixing length grows linearly in time. When the flow is channelling or in the single-finger exchange-flow regime (regimes II and IV) the mixing length tends to a constant value. As $\unicode[STIX]{x1D70E}$ is varied from $0$ (the homogeneous medium scenario), we see a transition from pure viscous fingering to channelling and fingering behaviour (figure 12 e). As $R$ is varied from $0$ (neutrally stable displacement scenario), we see a transition from pure channelling behaviour to both channelling and fingering (figure 12 f). Exactly which regimes occur depend on all four of the variables in this problem. In general we find fingering behaviour (regimes I and III) is more significant for larger values of $R$ and $Pe$ , while channelling is more significant for larger values of $\unicode[STIX]{x1D70E}$ . Also, fingering across layers (regimes III and IV) is only relevant when $n>1$ , and is most notable when the length scale of the permeability variations is small, $n\gg 1$ .

In the following sections we review the different regimes briefly and discuss some simple models for the evolution of the transversely averaged concentration. We note that the following is not an exhaustive description of all possible behaviour, and we leave an exhaustive description of the interplay between viscous fingering and permeability layering for future work.

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 1986), 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 1997b ; Shahnazari et al. 2018). 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. 2018, 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. (2018) 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. 2018). 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}$ .

7 Discussion and conclusions

In this paper, we have examined miscible displacements in heterogeneous porous media. The main goal throughout this work has been to better understand the structure and evolution of the concentration field during stable and unstable displacement processes through the use of high-resolution numerical simulations. Motivated by the fact that many geological formations consist of layered sedimentary sequences, we considered porous media with a permeability structure that varies perpendicular to the flow direction. We considered cases where the injected fluid is equally viscous, more viscous or less viscous than the ambient fluid. In general we find that the flow evolves through three main flow regimes. Figure 17 summarizes these different possible regimes in the cases of viscously dominated stable displacements (figure 17 a), permeability-dominated displacements (figure 17 b), and viscously dominated unstable displacements (figure 17 c). In each case the figure shows the instantaneous scaling exponent of the mixing length, $\unicode[STIX]{x1D6FF}$ , where $h=At^{\unicode[STIX]{x1D6FF}}$ , for different $Pe$ .

Figure 17. Representative plots of the scaling exponent of the mixing length, $\unicode[STIX]{x1D6FF}$ , found by locally fitting a power law of the form $h=At^{\unicode[STIX]{x1D6FF}}$ , for three different parameter sets: (a) $(R,\unicode[STIX]{x1D70E},n)=(-2,0.1,1)$ , demonstrating a viscously dominated stable displacement, as in § 5, (b) $(R,\unicode[STIX]{x1D70E},n)=(0.1,1,1)$ , demonstrating a permeability-dominated displacement, as in § 4 and (c) $(R,\unicode[STIX]{x1D70E},n)=(2,0.1,6)$ , demonstrating a viscously dominated unstable displacement, as in § 6. The regime boundaries (black lines) divide the early-time (regime 1), intermediate-time (regime 2) and late-time (regime 3) regimes. The dotted lines in (c) delineate the various intermediate-time regimes (§§ 6.1–6.3) and the dashed line delineates the hydrodynamically stable and unstable regions (Nijjer et al. 2018).

At early times (regime 1 in figure 17), the concentration field evolves through diffusion across an initially sharp interface ( $\unicode[STIX]{x1D6FF}=1/2$ ) and is independent of both the log-viscosity ratio, $R$ , and log-permeability ratio, $\unicode[STIX]{x1D70E}$ (see § 3.1). Once advection begins to outpace diffusion ( $t\sim O(1/Pe)$ ) the flow transitions to the intermediate-time regime. At intermediate times (regime 2 in figure 17), the interplay between viscosity and permeability variations leads to a range of possible dynamics ( $\unicode[STIX]{x1D6FF}\neq 1/2$ ). Finally, once the interface has become long and thin and transversely homogenized ( $t\sim O(Pe)$ ) the flow transitions to the late-time regime. At late times (regime 3 in figure 17), the flow becomes dominated by shear-enhanced dispersion and the concentration evolves diffusively again ( $\unicode[STIX]{x1D6FF}=1/2$ ).

At intermediate times (regime 2 in figure 17), the interplay between viscosity and permeability variations leads to different behaviour depending on the relative size of $R$ and $\unicode[STIX]{x1D70E}$ and whether the injected fluid is more viscous or less viscous than the ambient fluid. When permeability effects dominate ( $\unicode[STIX]{x1D70E}>|R|$ ; § 4); the viscosity modulates the effective permeability of the medium but otherwise evolves qualitatively in the same manner as when the two fluids have equal viscosity. When the injected fluid is more viscous than the ambient ( $|R|>\unicode[STIX]{x1D70E}$ and $R<0$ ; § 5), the viscosity contrast tends to prevent channelling at the interface and reduces spreading of the two fluids relative to the equal viscosity, $R=0$ case. When the injected fluid is less viscous than the ambient ( $|R|>\unicode[STIX]{x1D70E}$ and $R>0$ ; § 6), a number of different possible regimes were identified (regimes 2I, 2II, 2III and 2IV in figure 17 c). These different regimes arose due to the interplay between viscous fingering, which tends to cascade through a range of length scales, and the permeability structure, which has a fixed length scale. Overall, depending on which regime the flow is in, the permeability heterogeneity can either enhance or temper spreading relative to the uniform permeability case.

Figure 18. Long-time evolution of $\overline{c}(x,t)$ for $\unicode[STIX]{x1D70E}=0.3$ , $(Pe,n)=(100,1)$ and $t$ ranging from 10 to 1000. The viscosity ratio is (a) stabilizing $R=-1$ , (b) neutrally stable $R=0$ and (c) de-stabilizing $R=1$ . The asymptotic, viscosity-independent solution, equation (3.4) with (3.11) is given by the dashed line.

At late times (regime 3 in figure 17), the concentration evolves through shear-enhanced dispersion (§ 4.3), which asymptotically becomes independent of the viscosity ratio and only depends on the permeability structure. Figure 18 demonstrates this behaviour: it shows the evolution of the concentration field for a neutrally stable displacement, large stabilizing displacement and a large de-stabilizing displacement. Although spreading is initially hindered ( $R<0$ ), or enhanced ( $R>0$ ), at late times they both tend to the same viscosity-independent self-similar solution. This means that for processes that occur over very small length scales or very long time scales, the viscosity difference becomes insignificant and the permeability structure dictates the rate of spreading and mixing of the two fluids.

In summary, the flow evolves through three distinct regimes: an early-time regime dominated by longitudinal diffusion, an intermediate-time regime dominated by longitudinal advection and a late-time regime dominated by shear-enhanced dispersion. Informed by high-resolution numerical simulations, we have developed simple models that capture the dominant physics in each of the regimes, which provide an easy way of quantitatively predicting the average behaviour of these systems.