Stable and unstable miscible displacements in layered porous media

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.

Miscible displacements in layered porous media 469 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, 1988Zimmerman & Homsy 1991Jha, 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 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 µ 2 . Another fluid, with viscosity µ 1 , which is fully miscible with the first, is injected at a constant flux Q.
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.

Problem formulation
A schematic of the problem geometry is given in figure 1. We consider a semiinfinite, two-dimensional porous strip with finite width a. For simplicity, the porosity of the medium, φ, is constant, but the permeability, k = k(y), varies in the direction perpendicular to the flow. A fluid of viscosity µ 1 is injected into the medium at a constant volumetric flux Q, which is initially saturated with an ambient fluid of viscosity µ 2 . The two fluids are fully miscible.
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; Here 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, where R = log(µ 2 /µ 1 ).

Non-dimensionalization
We non-dimensionalize the governing equations by the width of the porous medium a, the injection flux Q, the mean permeabilityk, viscosity µ 2 , pressure µ 2 Q/k and Miscible displacements in layered porous media 471 convective time scale φa 2 /Q. We also change to a reference frame moving with the average speed of the injected fluid by settingx = x − Qt/a andũ = u − Q/a. Equations (2.1)-(2.4) become − (ũ + 1)µ = k ∂p ∂x , −vµ = k ∂p ∂y , (2.5a,b) (2.8) 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.

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 ln(k) = −σ cos(2πny) − ln (I 0 (σ )), (2.9) where 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 σ 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ŷ = ny,x = nx,t = n 2 t, in which case the flow evolves exactly as it would with n = 1, but with an effective Péclet numberPe = 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. 472 J. S. Nijjer, D. R. Hewitt and J. A. Neufeld

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) The upstream and downstream fluxes are zero (in the moving frame) and the transverse velocity vanishes in the far field so that, (2.14) We initialize the concentration field to have a step jump, 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 c(x, t) = 1 0 c dy and the mixing length h(t), defined below. To understand how 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 = c − c, ∂c ∂t and where f = 1 0 f dy 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 Miscible displacements in layered porous media 473 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) Note that we use this measure instead of the more common definition, h * = x| 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.

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) = (∂Ψ /∂y, −∂Ψ /∂x), to give The transverse velocity boundary conditions (2.11) become (2.20) We imposed the boundary conditions in the flow direction (2.12)-(2.14), at a finite distance x = ±Γ so that, where Γ was chosen to be sufficiently large such that the mixing region was always far from the boundaries. In fact, we increased Γ 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 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), u · ∇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.

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 2a). 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 Miscible displacements in layered porous media 475 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 In the case of sinusoidally varying log-permeability (2.9), the velocity is, Given this fixed, known velocity the concentration simply evolves via the advectiondiffusion 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 Using the initial and boundary conditions in § 2.3 the concentration evolves self-similarly as which holds at all times when the permeability is homogeneous (k = 1). The mixing length grows like h ∼ t 1/2 and can be calculated explicitly by substituting (3.4) into (2.18) (figure 3a). 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( ut) (where 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 u 2 ). 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 Plot of the transversely averaged concentration, c(x, t) versus x/t for (R, Pe, n) = (0, 500, 1), σ 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, u ∼ sinh(σ )/I 0 (σ ).

Intermediate-time behaviour: advection
given the initial condition (2.15). The transversely averaged concentration, c(x, t) can be calculated by averaging (3.6), The model gives good agreement with the numerical simulations (figure 3b) 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 ∼ ut (figures 3a and 4b).
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 c and a long, thin mixing zone, (2.17) reduces to while the transversely averaged concentration still evolves according to (2.16). Given that c is independent of y, we integrate this equation twice and impose periodicity and Substituting and solving for the convective flux in (2.16), using the expression for the velocity (3.1), leads to This convective flux can be written in the form of an effective diffusivity such that (2.16) reduces to is the shear-enhanced dispersivity, which only depends on the permeability structure (cf. Van den Broeck & Mazo 1983). For our choice of sinusoidally varying log-permeability, this integral cannot be solved analytically, but is instead integrated numerically for varying σ and plotted in figure 4(a). In appendix B, we derive the asymptotic limits of the shear-enhanced dispersivity for large and small σ (given as dashed and dot-dashed lines respectively in figure 4a). The solution to (3.11a) is again the similarity solution (3.4), but now with a modified Péclet number Pe * (3.11b) (see figure 4c). When the total dispersivity is dominated by the shear-enhanced dispersivity, Pe 2 S 1, the effective dispersion scales like Pe * ∼ 1/(Pe 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, intermediatetime advection and late-time shear-enhanced dispersion.

Small viscosity variations |R| < σ
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| < σ .
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, σ |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 σ = 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 5a,c,e); and at late times the flow is dominated by shear-enhanced dispersion (figure 5b,d,f ). The main difference between flows where R = 0 and R = 0 is that at intermediate times the interface is either Miscible displacements in layered porous media 479 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.

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 Combining (4.1) with the fact that the flux vanishes in any transverse slice in the moving frame, the velocity can be written as 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) = µ −1 / 1 0 µ −1 dy − 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 , 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 = c(x) + c (x, y), equation (4.2) becomes independent of the average concentration and only depends on the transverse variations, In the case of sinusoidal log-permeability variations, equation (4.3) is  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,ḣ = dh/dt, 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 6b), although it differs appreciably from the R = 0 case at the downstream tips (cf. snapshots in figure 5a,c,e).
In this regime, we first note that diffusion is negligible. This results in concentration deviations that are almost exactly either c = −c or c = 1 − c (figure 6c). 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 Miscible displacements in layered porous media 481 the domain aligned with the permeability structure and with magnitude 1/2 (dashed black line in figure 6c). Substituting these average deviations into (4.4), the mean streamwise velocity reduces to u = e −(σ +R/2) cos (2πy) 1 0 e −(σ +R/2) cos(2πs) ds − 1. (4.5) This approximate model results in intermediate-time dynamics equivalent to the neutrally stable case, but with an effective log-permeability ratio σ eff = σ + R/2. Given (4.5), we can calculate h as in § 3.2, and extractḣ by fitting a linear profile h =ḣt.
Modelling the effective permeability in this way gives reasonably good agreement with the numerical simulations (figure 6d), 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 = e σ +R /I 0 (σ ), and slower, u = e −(σ +R) /I 0 (σ ), 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 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,c = Rc , such that (4.6) becomes where β = R Pe ∂c/∂x = Pe ∂ ln µ(c)/∂x, incorporates all of the parameters in the problem, and can be thought of as a rescaled bulk concentration or viscosity gradient. Since β is independent of y, we can solve forc by integrating (4.7) twice. We perform this integration numerically, although the limits of small β, which is relevant here and is considered in § 4.3.1, and large β, which we find to be useful and is considered later in § 5.1, can be treated analytically. Having solved forc, we then calculate the shear-enhanced dispersivity,  (4.7). The leading-order small-|β| asymptotic behaviour, which is equal to when R = 0, is given by the dashed red line and the next-order corrections in β are given by the dot-dashed blue lines, equation (4.12).
(cf. (3.12)). Solutions of S for σ = 1 and R > 0 and R < 0 are given in figure 7. When R = 0, (4.7) reduces to (3.8) and S is exactly as described by (3.12), a permeabilitydependent 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 β ∝ ∂c/∂x, and diffusion always causes concentration gradients to diminish in time, β generally decreases over time. When |β| is large (i.e. the viscosity gradient is large, as at early times), and R > 0, S diverges, whereas when R < 0, 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 |β| is small (i.e. the viscosity gradient is small, as at late times), S becomes independent of β 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. When |R| < σ , we need to only consider the small-|β| limit because by the time the flow reaches the late-time regime, β 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 u), or equivalently, the concentration gradient is ∂c/∂x = O(1/Pe u). This means that when the flow transitions to the late-time regime, |β| = O(R/ u) < O(1).
(4.9a,b) Given that the concentration deviations must satisfy periodicity and have vanishing mean, we find thatc 1 = I(k − 1) − I(k − 1), (4.10) c 2 = I(kc 1 ) − I(kc 1 ) − kc 1 I(k) + kc 1 I(k), (4.11) where I( f ) = y 0 ζ 0 f dη dζ for a given function f and, as before, the overbar refers to a transverse average. The shear-enhanced dispersivity, S, is thus (4.12) The leading-order contribution to 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 dotdashed blue lines.

Comparison to numerical simulations
We now use the calculated shear-enhanced dispersivity to determine how c evolves in time. Allowing S to vary in space, (2.16) yields a nonlinear diffusion equation for c, where β(x) = R Pe ∂c/∂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 S (by solving (4.7) and calculating (4.8)), and the small-|β| 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(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 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).

Large stabilizing viscosity variations |R| > σ , 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(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 9a). This is because, for large |∂c/∂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 9f ). 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 9b,c). Eventually the interface becomes very diffuse and the velocity becomes predominantly longitudinal (figure 9d) as in the cases considered in previous sections. In contrast to those cases, however, where the concentration deviations, (c ), 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.

Concentration model
As found in § 4.3, the shear-enhanced dispersivity S is given by (4.8) where β = R Pe ∂c/∂x andc 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 1). In particular, |β| is not just small but can take on any value. We thus first consider the large-β limit.
We start by expandingc in powers of 1/β,c =c 0 +c 1 /β + O(1/β 2 ). Substituting into (4.7) expanding and equating different powers of β, we find where we have again used the fact that 1 0c 0 = 1 0c 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 In the case of sinusoidally varying log-permeability, this limit corresponds to 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, where S * = 1 0 [ y 0 (k − 1) dη] 2 dy 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). for (R, σ , 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 S (found by solving (4.7)), or (5.6), are given by the dashed and dot-dashed black lines respectively.

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 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 σ in figure 11(a). The very good agreement between the reduced-order model and the full simulations over a range of σ 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. 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 > σ . In the absence of any permeability variations, σ = 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.

Miscible displacements in layered porous media
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 12a). In the second 488 J. S. Nijjer, D. R. Hewitt and J. A. Neufeld regime (II), the flow follows the imposed layered structure while diffusing across the layers causing the flow to slow down (figure 12b). 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 12c). These fingers then also coarsen, leading to a fourth regime (IV) where a single pair of counter-propagating fingers remain (figure 12d). 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 12e,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 σ is varied from 0 (the homogeneous medium scenario), we see a transition from pure viscous fingering to channelling and fingering behaviour (figure 12e). As R is varied from 0 (neutrally stable displacement scenario), we see a transition from pure channelling behaviour to both channelling and fingering (figure 12f ). 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 σ . 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 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.

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 ∼ 1/RPe (Tan & Homsy 1986), is small compared to the length scale imposed by the permeability y ∼ 1/n. Hence the flow always fingers for sufficiently large Pe and is also observed to be enhanced when σ 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 < σ , so long as Pe is sufficiently large. Figure 13(a,b) shows snapshots of the concentration field for σ = 0.1 and σ = 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 13c). 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 13d), and 1 0 c 2 dy versus the similarity variable x/t for σ 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. 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 σ for σ 0.1 ( figure 13c,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 σ 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.

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 12b). Figure 14(a-d) shows a sequence of snapshots of the concentration field with overlain streamlines for (R, σ , 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 14f ). As time progresses, and the fluids become more mixed, the effect of the viscosity reduces and the This behaviour resembles the late-time regime in the miscible viscous-fingering instability (i.e. when σ = 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 15a). 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, α, to the results of the 2-D numerical simulations (figure 15b). As the viscosity ratio and permeability variance are increased α 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 16a,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 (σ = 0).
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 16c). 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 16d). 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 < σ . This can occur when the system reaches the late-time regime and the interface is still hydrodynamically unstable. For the interface  (Nijjer et al. 2018).
to be stable, the width of the mixing zone must be h ∼ RPe (Nijjer et al. 2018). The mixing zone at the start of the late-time regime in the permeability-dominated scenario is h ∼ Pe u/n 2 ∼ Peσ /n 2 and so the late-time state may be unstable to further coarsening via a viscous instability if R > σ 2 /n 2 .

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 17a), permeability-dominated displacements (figure 17b), and viscously dominated unstable displacements (figure 17c). In each case the figure shows the instantaneous scaling exponent of the mixing length, δ, where h = At δ , for different Pe.
At early times (regime 1 in figure 17), the concentration field evolves through diffusion across an initially sharp interface (δ = 1/2) and is independent of both the log-viscosity ratio, R, and log-permeability ratio, σ (see § 3.1). Once advection begins to outpace diffusion of R and σ and whether the injected fluid is more viscous or less viscous than the ambient fluid. When permeability effects dominate (σ > |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| > σ 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| > σ and R > 0; § 6), a number of different possible regimes were identified (regimes 2I, 2II, 2III and 2IV in figure 17c). 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. At late times (regime 3 in figure 17), the concentration evolves through shearenhanced 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. It can be seen that h calculated using (A 1) shows no discernible difference between the two simulations, whereas h calculated using (2.18) is noticeably different at intermediate times. This is because the simulations shown in figure 19 correspond to large stabilizing viscosities whose concentration profiles are composed of a sharp gradient at the interface, which depends on the viscosity ratio, and long tails upstream and downstream, which are independent of the viscosity ratio. Since (A 1) only picks up the growth of the tail regions, it is independent of R, whereas (2.18) accounts for the evolution of the sharp concentration gradient and hence depends on R. We therefore use (2.18) to measure the mixing length due to this sensitivity to the structure of the concentration field. When the variations in the permeability are large (σ 1), the fluid velocity diverges at y = 1/2 and tends to −1 everywhere else. The charactersitic spreading velocity is Appendix C. Concentration models for R > σ C.1. Regime I: fingering within layers To model the nonlinear evolution of the fingered region, we again assume the flow is in vertical flow equilibrium and that the streamwise velocity is given by (4.3). However, unlike in § § 3.2 and 4.2 where the velocity depends on the permeability, here we assume the velocity only depends on the viscosity, Next we consider η f (x) forward-propagating fingers of width w f (x) and η b (x) backward-propagating fingers of width w b (x). We make the simplifying assumption that the forward-and backward-propagating fingers have uniform concentration c = 1 and c = 0 respectively, but allow the viscosity to vary inside the finger. We have the additional constraints that the transversely averaged concentration is equal to the proportion of forward-propagating fingers η f w f = c and the total area of the fingers adds up to one, η f w f + η b w b = 1.
In a homogeneous medium (σ = 0), the forward-and backward-propagating fingers have viscosities µ(y) ≈ e R(1−y 2 /w 2 f ) and µ(y) ≈ e R(y 2 /w 2 f ) across the fingers respectively (see Nijjer et al. 2018). When σ is non-zero, the backward-propagating fingers are broad and coincide with the low permeability layers. In this case, the viscosity in the fingers is nearly uniform, µ(y) ≈ 1.
Substituting our assumptions for the concentration and viscosity into (C 2) gives, where The transversely averaged concentration evolves in exactly the same manner as in the homogeneous case but with a larger effective viscosity (cf. Nijjer et al. 2018). This model prediction is given by the dot-dashed line in figure 13(c) and gives good quantitative agreement with the numerical simulations when σ 0.1.
C.2. Regime II: stable channelling This regime is characterized by a linear background gradient with superimposed deviations that decay. In § 5 we found that in the limit of large aspect ratio and small, quasi-statically evolving deviations, (2.17) reduced to (4.6). Here, in the channelling regime, the deviations decay exponentially, independently of the transversely averaged concentration. We therefore generalize the analysis of § 4.3 by allowing the deviations to evolve in time. Equation (2.16) remains the same, but (2.17) becomes ∂c ∂t + u ∂c ∂x = 1 Pe ∂ 2 c ∂y 2 . (C 6) Again, the fluid flow is assumed to be in vertical flow equilibrium and so the velocity is given by (4.4). For simplicity, we will also assume σ 1, although the method below can be generalized to σ = O(1). Since c c = O(1), (4.4) can be written, to leading order, as u = Rc − σ cos(2nπy), (C 7) and (C 6) becomes ∂c ∂t + Rc ∂c ∂x − 1 Pe ∂ 2 c ∂y 2 = σ cos(2nπy) ∂c ∂x .
(C 8) When σ = 0 this is exactly the late-time, single-finger exchange-flow behaviour of the miscible viscous-fingering instability (Nijjer et al. 2018). Solving (C 8) using separation of variables yields c = −σ ∂c/∂x γ + Ke −γ t cos(2nπy), (C 9) where γ = R∂c/∂x + 4n 2 π 2 /Pe. The constant of integration K corresponds to the initial conditions at the onset of this regime and incorporates the nonlinear earlytime spreading. A simple approximation would be that K ≈ 1/2 corresponding to the maximum amplitude of sinusoidally varying c . Depending on which term is dominant in (C 9), one gets one of two limiting solutions to (2.16). When t is small (intermediate times), the second term in (C 9) dominates over the first. The deviations are longitudinally uniform and decay exponentially, and the solution of (2.16) for c is c = 1 2 − αx.
(C 10) This tendency to a steady linear profile can be seen in figure 15(a). When t is large (late times), the first term dominates over the second. In this case the deviations decay quasi-steadily with the transversely averaged concentration. This case is exactly the small-β limit discussed in § 4.3.1 in the limit of small σ , and is dominated by shearenhanced dispersion. In this case c tends to the self-similar error-function solution (3.4) with effective diffusivity (3.11). In addition, if σ > R, the first term, which scales like O(σ /R) (given ∂c/∂x ∼ 1/RPe), is always larger than Ke −γ t , consistent with the fact that this channelling behaviour is only present when viscosity variations are more important than permeability variations, R > σ . Thus far the model is unable to predict the slope of the well-mixed interior region, α. To close the model, we relate the convective flux to the mixing that occurs. We consider a control volume containing everything to the right of the midplane. Neglecting longitudinal diffusion, the total convective flux into this control volume must balance the net change in concentration where T is a timescale marking the end of this regime. Assuming the flow is always in this channelling regime, we can substitute (C 7) and (C 8) into (C 11) to attain an implicit equation for α. This control volume approach inherently depends on the choice of T; as T → ∞, the transversely averaged concentration evolves via shearenhanced dispersion and the estimate for α → 0. We therefore select T at the boundary between the intermediate-and late-time regimes such that the concentration profile is linearly well mixed in the interior and has not transitioned to the late-time diffusive solution. For our simulations, we find the choice T = 200 ≈ 3/γ appropriate, although we note that predictions are broadly insensitive to the choice of T as long as T ∼ O(1/γ ).