Stability Results on Radial Porous Media and Hele-Shaw Flows with Variable Viscosity Between Two Moving Interfaces

TL;DR

This paper conducts a linear stability analysis of three-layer radial flows with variable viscosity, deriving bounds on eigenvalues and identifying optimal viscosity profiles to understand flow stability.

Contribution

It introduces a novel eigenvalue problem with time-dependent coefficients for analyzing flow stability with variable viscosity in multi-layer radial flows.

Findings

Derived upper bounds on the eigenvalue spectrum.

Characterized eigenvalues and eigenfunctions in $L^2$ space.

Numerically computed eigenvalues and identified optimal viscosity profiles.

Abstract

We perform a linear stability analysis of three-layer radial porous media and Hele-Shaw flows with variable viscosity in the middle layer. A nonlinear change of variables results in an eigenvalue problem that has time-dependent coefficients and eigenvalue-dependent boundary conditions. We study this eigenvalue problem and find upper bounds on the spectrum. We also give a characterization of the eigenvalues and prescribe a measure for which the eigenfunctions are complete in the corresponding $L^2$ space. The limit as the viscous gradient goes to zero is compared with previous results on multi-layer radial flows. We then numerically compute the eigenvalues and obtain, among other results, optimal profiles within certain classes of functions.

Figures

Equations

$$\nabla \cdot \mathbf{u}=0,\nabla p=-\mu \mathbf{u},\frac{\partial \mu }{\partial t}+\mathbf{u}\cdot \nabla \mu =0,\text{for }r\ne 0.$$

$$\zeta =\frac{r^2-R_0^2(t)}{R_2^2(t)-R_0^2(t)}=\frac{r^2-R_0^2(t)}{R_2^2(0)},$$

$$\alpha =\theta ,$$

$$\tau =t.$$

$$\left.\begin{array}[]{l}\frac{\partial\tilde{u}_{\zeta}}{\partial\zeta}+\frac{1}{\zeta}\frac{\partial\tilde{u}_{\alpha}}{\partial\alpha}=0\\ \frac{\partial\tilde{p}}{\partial\zeta}=-\frac{R_{2}^{4}(0)}{4(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau))}\mu\tilde{u}_{\zeta}-\frac{QR_{2}^{2}(0)}{4\pi(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau))}\tilde{\mu}\\ \frac{\partial\tilde{p}}{\partial\alpha}=-\frac{\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau)}{\zeta}\mu\tilde{u}_{\alpha}\\ \frac{\partial\tilde{\mu}}{\partial\tau}+\tilde{u}_{\zeta}\frac{\partial\mu}{\partial\zeta}=0.\end{array}\right\}$$

$$\left(\tilde{u}_{\zeta},\tilde{u}_{\alpha},\tilde{p},\tilde{\mu}\right)=\big{(}f(\zeta,\tau),\delta(\zeta,\tau),\psi(\zeta,\tau),\phi(\zeta,\tau)\big{)}e^{in\alpha}.$$

$$\frac{\partial \varphi (\zeta ,\tau )}{\partial \tau }=-\frac{d\mu }{d\zeta }f(\zeta ,\tau ),$$

$$\begin{array}{rl}& \frac{\partial }{\partial \zeta }\left\{\left(\zeta R_2^2(0)+R_0^2(\tau )\right)\mu \frac{\partial f(\zeta ,\tau )}{\partial \zeta }\right\}-\frac{n^2{R}_2^4(0)}{4(\zeta R_2^2(0)+R_0^2(\tau ))}\mu f(\zeta ,\tau )\\ =& \frac{Q{n}^2{R}_2^2(0)}{4\pi (\zeta R_2^2(0)+R_0^2(\tau ))}\varphi (\zeta ,\tau ).\end{array}$$

$$f(\zeta ,\tau )=\widetilde{C_1}{\left(\zeta R_2^2(0)+R_0^2(\tau )\right)}^{\frac{n}{2}}+\widetilde{C_2}{\left(\zeta R_2^2(0)+R_0^2(\tau )\right)}^{-\frac{n}{2}}.$$

$$C_n^{\prime }(\tau )=f({\zeta }_1,\tau ).$$

$$D_n^{\prime }(\tau )=f(1,\tau ).$$

$$\frac{2R_{1}^{2}(\tau)}{R_{2}^{2}(0)}\Big{(}\tilde{p}^{+}(\zeta_{1})-\tilde{p}^{-}(\zeta_{1})\Big{)}-\eta_{i}\frac{Q}{2\pi}\Big{(}\mu(\zeta_{1})-\mu_{i}\Big{)}=T_{1}\left(\frac{\eta_{i}+\frac{\partial^{2}\eta_{i}}{\partial\alpha^{2}}}{R_{1}(\tau)}\right),$$

$$\frac{2R_{1}^{2}(\tau)}{R_{2}^{2}(0)}\Big{(}\mu_{i}(f^{-})^{\prime}(\zeta_{1},\tau)-\mu(\zeta_{1})(f^{+})^{\prime}(\zeta_{1},\tau)\Big{)}=\left\{\frac{Qn^{2}}{2\pi R_{1}^{2}(\tau)}\Big{(}\mu(\zeta_{1})-\mu_{i}\Big{)}-T_{1}\frac{n^{4}-n^{2}}{R_{1}^{3}(\tau)}\right\}C_{n}(\tau).$$

$$(f^{-}{)}^{\prime }({\zeta }_1,\tau )=\frac{n{R}_2^2(0)}{2R_1^2(\tau )}f({\zeta }_1,\tau ).$$

$$n\mu_{i}f(\zeta_{1},\tau)-\frac{2R_{1}^{2}(\tau)}{R_{2}^{2}(0)}\mu(\zeta_{1})(f^{+})^{\prime}(\zeta_{1},\tau)=\left\{\frac{Qn^{2}}{2\pi R_{1}^{2}(\tau)}\Big{(}\mu(\zeta_{1})-\mu_{i}\Big{)}-T_{1}\frac{n^{4}-n^{2}}{R_{1}^{3}(\tau)}\right\}C_{n}(\tau).$$

$$C_n^{\prime }(\tau )=\frac{2R_1^2(\tau )}{n{R}_2^2(0)}\frac{\mu ({\zeta }_1)}{{\mu }_i}(f^{+}{)}^{\prime }({\zeta }_1,\tau )+\frac{F_1}{{\mu }_i}{C}_n(\tau ),$$

$$F_{1}=\frac{Qn}{2\pi R_{1}^{2}(\tau)}\Big{(}\mu(\zeta_{1})-\mu_{i}\Big{)}-T_{1}\frac{n^{3}-n}{R_{1}^{3}(\tau)}.$$

$$D_n^{\prime }(\tau )=-\frac{2R_2^2(\tau )}{n{R}_2^2(0)}\frac{\mu (1)}{{\mu }_o}(f^{-}{)}^{\prime }(1,\tau )+\frac{F_2}{{\mu }_o}{D}_n(\tau ),$$

$$F_{2}=\frac{Qn}{2\pi R_{2}^{2}(\tau)}\Big{(}\mu_{o}-\mu(1)\Big{)}-T_{2}\frac{n^{3}-n}{R_{2}^{3}(\tau)},$$

$$\left.\begin{array}[]{l l}\frac{\partial\phi(\zeta,\tau)}{\partial\tau}=-\frac{d\mu}{d\zeta}f(\zeta,\tau),\\ \frac{\partial}{\partial\zeta}\left\{\left(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau)\right)\mu\frac{\partial f(\zeta,\tau)}{\partial\zeta}\right\}-\frac{n^{2}R_{2}^{4}(0)}{4(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau))}\mu f(\zeta,\tau)=\frac{Qn^{2}R_{2}^{2}(0)\phi(\zeta,\tau)}{4\pi(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau))},&\\ C_{n}^{\prime}(\tau)=\frac{2R_{1}^{2}(\tau)}{nR_{2}^{2}(0)}\frac{\mu(\zeta_{1})}{\mu_{i}}f^{\prime}(\zeta_{1})+\frac{F_{1}}{\mu_{i}}C_{n}(\tau),&\\ D_{n}^{\prime}(\tau)=-\frac{2R_{2}^{2}(\tau)}{nR_{2}^{2}(0)}\frac{\mu(1)}{\mu_{o}}f^{\prime}(1)+\frac{F_{2}}{\mu_{o}}D_{n}(\tau),&\end{array}\right\}$$

$$\left.\begin{array}[]{l}\frac{\partial\phi(\zeta,\tau)}{\partial\tau}=\sigma(\tau)\phi(\zeta,\tau),\\ C_{n}^{\prime}(\tau)=\sigma(\tau)C_{n}(\tau),\\ D_{n}^{\prime}(\tau)=\sigma(\tau)D_{n}(\tau),\end{array}\right\}$$

$$\left.\begin{array}[]{l l}\Big{(}\left(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau)\right)\mu f^{\prime}(\zeta)\Big{)}^{\prime}-\frac{n^{2}R_{2}^{4}(0)}{4(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau))}\mu f(\zeta)=-\frac{Qn^{2}R_{2}^{2}(0)}{4\pi(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau))}\frac{1}{\sigma(\tau)}\frac{d\mu}{d\zeta}f(\zeta),\\ \frac{2R_{1}^{2}(\tau)}{nR_{2}^{2}(0)}\mu(\zeta_{1})f^{\prime}(\zeta_{1})=\left(\mu_{i}-\frac{F_{1}}{\sigma(\tau)}\right)f(\zeta_{1}),&\\ -\frac{2R_{2}^{2}(\tau)}{nR_{2}^{2}(0)}\mu(1)f^{\prime}(1)=\left(\mu_{o}-\frac{F_{2}}{\sigma(\tau)}\right)f(1).&\end{array}\right\}$$

$$r=\sqrt{\zeta R_2^2(0)+R_0^2(t)}=\sqrt{({\zeta }_1+C_n(\tau )e^{in\alpha })R_2^2(0)+R_0^2(t)}.$$

$$r=R_1(\tau )+\frac{R_2^2(0)}{2R_1(\tau )}{C}_n(\tau )e^{in\alpha }+\mathcal{O}(C_n^2(\tau )).$$

$$A_n(t)=\frac{R_2^2(0)}{2R_1(\tau )}{C}_n(\tau ).$$

$$\frac{A_n^{\prime }(t)}{A_n(t)}=\frac{C_n^{\prime }(\tau )}{C_n(\tau )}-\frac{Q}{2\pi R_1^2(t)}.$$

$$\frac{B_n^{\prime }(t)}{B_n(t)}=\frac{D_n^{\prime }(\tau )}{D_n(\tau )}-\frac{Q}{2\pi R_2^2(t)},$$

$$\zeta =\frac{r^2-R_0^2(t)}{R^2(0)}.$$

$$D_n^{\prime }(\tau )=f(1,\tau ),$$

$$\frac{2R^{2}(\tau)}{R^{2}(0)}\Big{(}-\mu_{o}(f^{+})^{\prime}(1,\tau)+\mu_{i}(f^{-})^{\prime}(1,\tau)\Big{)}=\left\{\frac{Qn^{2}}{2\pi R^{2}(\tau)}\Big{(}\mu_{o}-\mu_{i}\Big{)}-T\frac{n^{4}-n^{2}}{R^{3}(\tau)}\right\}D_{n}(\tau),$$

$$f(\zeta,\tau)=\widetilde{C_{1}}\Big{(}\zeta R^{2}(0)+R_{0}^{2}(\tau)\Big{)}^{\frac{n}{2}},\qquad\zeta<1,$$

$$f(\zeta,\tau)=\widetilde{C_{2}}\Big{(}\zeta R^{2}(0)+R_{0}^{2}(\tau)\Big{)}^{-\frac{n}{2}},\qquad\zeta>1.$$

Full text

1 Introduction

Saffman-Taylor instability is an instability that occurs when a less viscous fluid drives a more viscous fluid in a porous medium. This phenomenon occurs in many applications including carbon sequestration, filtration, hydrology, and petrology. One important application is oil recovery in which it is a common practice to use water to displace oil. Oil recovery was the application driving the seminal work on this type of instability by Saffman & Taylor, (1958). In order to simplify their experiments, Saffman and Taylor studied the instability in the context of Hele-Shaw flows - flows between two parallel plates with a small gap between them. Hele-Shaw flows are a good model of porous media flows because they are also governed by Darcy’s Law. It is now well known that a positive interfacial viscous jump in the direction of rectilinear flow produces an unstable flow with interfacial tension stabilizing short waves.

In applications in which it is advantageous to suppress the instability, one strategy is to use intermediate fluids which are more viscous than the displacing fluid but less viscous than the displaced fluid in order to limit the size of the viscous jumps. In the case that the interfacial tension at the interfaces is comparable to that of two-layer flow, this strategy is effective. This has been demonstrated by Daripa, (2008) in the case of rectilinear flow using upper bounds on the growth rates which were also used to give stabilization criteria.

Another approach to limiting the instability is using a variable viscosity fluid as the intermediate fluid. This can be achieved, for example, in chemical enhanced oil recovery in which polymer is used to increase the viscosity of the displacing fluid and variable polymer concentration leads to variable viscosity. A viscous gradient allows for even smaller viscous jumps at the interfaces than a constant viscosity fluid but the layer itself becomes unstable. There are several studies on the stability of multi-layer variable viscosity porous media and Hele-Shaw flows in a rectilinear geometry. Gorell & Homsy, (1983) first theoretically studied the stability of three-layer flows in such a geometry in which the middle fluid has variable viscosity. However, they studied the restricted case in which the trailing interface is a miscible interface with no interfacial tension. Daripa & Pasa, (2006) dropped this restriction and studied the case of three-layer variable viscosity flow with both interfaces immiscible.

Saffman-Taylor instability is also studied in a radial flow geometry in which there is a point source in the center of the flow and the fluid moves outward radially with circular interfaces between fluids. Radial flow is one of several cases investigated by Muskat, (1934, 1946) for the two-layer problem but in the case of zero interfacial tension. Paterson, (1981) later performed a linear stability analysis for two-layer radial flow with interfacial tension. There are relatively few works on the stability of multi-layer radial flows. Cardoso & Woods, (1995) studied the stability of three-layer radial flows in the limiting case in which the inner interface is completely stable and looked at the break up of the middle layer into drops. Beeson-Jones & Woods, (2015) analyzed general three-layer flows, and Gin & Daripa, (2015) performed a linear stability analysis for an arbitrary number of fluid layers in radial geometry.

To date, there are no known studies of multi-layer radial flows with variable viscosity. The development of the theory for radial flow lags behind that of rectilinear flow because of the challenges due to the time-dependence of the problem. In particular, the curvature of the interfaces, the length of the middle layer(s) of fluid, and the spatially dependent viscous profile are all time-dependent. Previous stability studies of radial Hele-Shaw flows (Anjos et al.,, 2015; Beeson-Jones & Woods,, 2015, 2017; Cardoso & Woods,, 1995; Dallaston & McCue,, 2013; Kim et al.,, 2009; Paterson,, 1981) allow for a time-dependent growth rate of instabilities, typically making use of a quasi-steady-state approximation. Despite this approximation, the stability results have often shown agreement with both physical and numerical experiments. In what follows, we keep with this tradition and use a quasi-steady-state approximation to derive an eigenvalue problem that has time-dependent coefficients and thus time-dependent eigenvalues. Additionally, the eigenvalue problem has eigenvalue-dependent boundary conditions which make the study of the eigenvalue problem mathematically interesting and challenging. The eigenvalues depend on many different parameters including the viscous profiles of the fluid layers, the interfacial tension at each interface, the curvature of the interfaces, and the fluid injection rate. Therefore, a numerical exploration of the vast parameter space is necessary. It is also important for applications in chemical EOR to investigate which viscous profiles minimize the instability for a given set of parameters, as is done in by Daripa & Ding, (2012) in the rectilinear case.

It is worth pointing out another significant motivation behind this study. The simplest model to study incompressible porous media flow is to use the Hele-Shaw model which is based on the incompressibility condition and the Darcy’s law in each of the phases. This model does not allow any mixing between the oil and water phases macroscopically and maintains sharp interfaces between the phases. Another model, called the Buckley-Leverett model, builds on this Hele-Shaw model by adding a saturation equation, a nonlinear hyperbolic conservation law with non-convex flux function which allows mixing between the phases (macroscopically) due to rarefaction waves behind the leading saturation front sweeping the oil ahead (Daripa et al.,, 1988). The rarefaction waves creates a viscous profile behind the front with viscosity gradually increasing towards the moving front. This mixing region with a viscous profile is usually finite in length which grows with time. The study of the stability of such composite solutions to such conservation laws is relevant for porous media but much too difficult. An appropriate model for this problem is the one under study which models the effect of rarefaction waves with a viscous profile between two interfaces and the shock front with a material interface having appropriate viscous jump at the interface. Such a study was partially carried out by Daripa, (2008) in the rectilinear geometry but the present study in the radial geometry is much harder as mentioned above and is also more relevant for porous media flow.

The paper is laid out as follows. In section 2 we perform a linear stability analysis of a point source driven three-layer radial Hele-Shaw flow in which fluid between the two interfaces has a variable viscous profile. We use a time-dependent coordinate transformation to freeze the basic motion of the two interfaces in this new coordinate system and derive the associated eigenvalue problem in this new coordinate system. The growth rate of disturbances in the transformed coordinate system is related to the physical growth of disturbances of the interfaces in section 3. Section 4 gives the restriction of the problem to the case of constant viscosity. In section 5, upper bounds on the growth rate are derived using the variational form of the problem. The nature of the eigenvalues and the completeness of the eigenfunctions are investigated in section 6. Numerical evaluation of the eigenvalues and their dependence on certain physical parameters are given in section 7, and then we conclude in section 8.

2 Preliminaries

We consider a radial Hele-Shaw flow consisting of three regions of incompressible, immiscible fluid. By averaging across the gap, we may consider a two-dimensional flow domain in polar coordinates, $\Omega:=(r,\theta)=\mathbb{R}^{2}$. The least viscous fluid with constant viscosity $\mu_{i}$ is injected into the center of the cell at a constant injection rate, $Q$. The most viscous fluid, with constant viscosity $\mu_{o}$, is the outermost fluid. The middle layer fluid has a smooth, axisymmetric viscous profile $\mu(r)$ where $\mu_{i}<\mu(r)<\mu_{o}$. The fluid flow is governed by the following equations

[TABLE]

The first equation $\eqref{VariableViscosity:eq:main}_{1}$ is the continuity equation for incompressible flow, the second equation $\eqref{VariableViscosity:eq:main}_{2}$ is Darcy’s Law, and the third equation $\eqref{VariableViscosity:eq:main}_{3}$ is an advection equation for viscosity. Note that Darcy’s Law for Hele-Shaw flows contains a permeability term $K=b^{2}/12$, but for the sake of simplicity here we have scaled viscosity by this term. Therefore, in what follows $\mu$ denotes the modified viscosity. We start with the fluids separated by circular interfaces with radii $R_{1}(0)$ and $R_{2}(0)$, where $R_{1}(t)$ and $R_{2}(t)$ are the positions of the interfaces at time $t$. This set-up is shown in Figure 1.

The equations admit a simple basic solution in which all of the fluid moves outward radially with velocity $\mathbf{u}:=\left(u_{r},u_{\theta}\right)=\left(Q/(2\pi r),0\right)$. The interfaces remain circular and their radii are given by $R_{1}(t)=\sqrt{Qt/\pi+R_{1}(0)^{2}}$ and $R_{2}(t)=\sqrt{Qt/\pi+R_{2}(0)^{2}}$. The pressure, $p_{b}=p_{b}(r)$, may be obtained by integrating equation $\eqref{VariableViscosity:eq:main}_{2}$.

We define the quantity $R_{0}(t)=\sqrt{Qt/\pi}$ and define the following coordinate transformation:

[TABLE]

The basic solution in these coordinates is $(u_{\zeta},u_{\alpha})=(0,0)$ with the interfaces stationary at $\zeta=\zeta_{1}:=R_{1}^{2}(0)/R_{2}^{2}(0)$ and $\zeta=1$. $\mu=\mu(\zeta)$ is now independent of time.

We perturb this basic solution $\left(u_{\zeta}=0,u_{\alpha}=0,p_{b},\mu\right)$ by $\left(\tilde{u}_{\zeta},\tilde{u}_{\alpha},\tilde{p},\tilde{\mu}\right)$ where the disturbances are assumed to be small. The linearized equations which govern these disturbances are

[TABLE]

We use separation of variables and decompose the disturbances into Fourier modes in the $\alpha$ coordinate so the disturbances are of the form

[TABLE]

Using this ansatz in equation (5) yields the following relations:

[TABLE]

In the innermost and outermost layers, the viscosity is constant and therefore $\phi(\zeta,\tau)\equiv 0$. In these regions, the solution of (8) is of the form

[TABLE]

2.1 Interface Conditions

Recall that the inner interface in the $\zeta$-coordinate system is located at $\zeta=\zeta_{1}:=R_{1}^{2}(0)/R_{2}^{2}(0)$. Let the disturbance of this interface be of the form $\eta_{i}=C_{n}(\tau)e^{in\alpha}$. The linearized kinematic condition at the inner interface is given by

[TABLE]

The outer interface is located at $\zeta=1$. If the disturbance of the interface is of the form $\eta_{o}=D_{n}(\tau)e^{in\alpha}$, then the linearized kinematic condition is

[TABLE]

The linearized dynamic interface condition at the inner interface is given by

[TABLE]

where $T_{1}$ is the interfacial tension. Using the ansatz (6) and the system (5),

[TABLE]

When $\zeta<\zeta_{1}$, $f$ is of the form given by (9). When $\tau=0$, in order to avoid a singularity when $\zeta\to 0$, $f$ must be of the form $f(\zeta,\tau)=\widetilde{C_{1}}(\zeta R_{2}^{2}(0)+R_{0}^{2}(\tau))^{\frac{n}{2}}.$ We assume this also to be true for $\tau>0$. Then

[TABLE]

Using (13) in (12),

[TABLE]

Combining this with the kinematic interface condition (10),

[TABLE]

where $F_{1}$ is given by

[TABLE]

A similar procedure for the outer interface yields the interface condition

[TABLE]

where $F_{2}$ is given by

[TABLE]

and $T_{2}$ is the interfacial tension at the outer interface.

2.2 Eigenvalue problem

To this point we have the following system of equations where the field equations hold in the domain $(\zeta_{1},1)$:

[TABLE]

where we have dropped the superscripts “+” and “-”. Using a quasi-steady-state approximation (QSSA) in which $\tau$ (and hence $R(\tau)$) are frozen so that the functions $\phi(x,\tau)$, $C_{n}(\tau)$, and $D_{n}(\tau)$ experience short-time exponential growth satisfying

[TABLE]

for some growth rate $\sigma(\tau)$. Plugging (20) into (19) and using (10) and (11), $(f,\sigma)$ is a solution to the following eigenvalue problem in the domain $(\zeta_{1},1)$:

[TABLE]

The eigenvalues of system (21) are the time-dependent growth rates of the disturbances of the system. The QSSA allows for the considerable analysis and computation of growth rates that follows.

3 Relating the growth of interfacial disturbances in the $\zeta$-coordinates with the physical coordinates

We now relate the growth of the interfacial disturbances in the $\zeta$-coordinates to the same in the physical coordinate system. We start with the inner interface. Recall that in the transformed coordinates, the inner interface was disturbed by $C_{n}(\tau)e^{in\alpha}$. Therefore, it is located at $\zeta=\zeta_{1}+C_{n}(\tau)e^{in\alpha}$. Thus, the position of the interface in the physical coordinates is

[TABLE]

Expanding about $\zeta=\zeta_{1}$,

[TABLE]

If we write the disturbance in the physical coordinates as $A_{n}(t)e^{in\theta}$ (that is, the interface is located at $r=R_{1}(t)+A_{n}(t)e^{in\theta}$), then, within linear approximation,

[TABLE]

This implies that

[TABLE]

Following the same process, the growth of the outer interface is

[TABLE]

where the outer interface is located at $r=R_{2}(t)+B_{n}(t)e^{in\theta}$.

4 Constant Viscosity Fluids

We now consider the case in which all of the fluids have constant viscosity, first for two-layer flows and then for three-layer. We do this to demonstrate that the variable viscosity formulation can recover previous results in the constant viscosity limit. Through this process, we also present some new results on three-layer constant viscosity flows.

When there are only two fluids (i.e. one interface located at $r=R(t)$), the above analysis holds with the coordinate transformation

[TABLE]

In the new coordinates, the basic solution has the interface fixed at $\zeta=1$. Let $\mu_{i}$ denote the viscosity of the inner fluid and $\mu_{o}$ denote the viscosity of the outer fluid. Analogous to equations (11) and (12), the kinematic interface condition is

[TABLE]

and the dynamic interface condition is

[TABLE]

where $T$ is the interfacial tension and $D_{n}(\tau)$ is the amplitude of the disturbance of the interface with wave number $n$. Also, as stated in the derivation of the interface conditions above,

[TABLE]

and

[TABLE]

Using these in equations (25) and (26) gives the two-layer growth rate

[TABLE]

This is an expression for the growth rate of the disturbance of the interface in the $\zeta$-coordinate system. This problem can be solved in the original $r$-coordinate system, and the result is a classic one (Paterson,, 1981). We recall this result, which has been reproduced using our current notation by Gin & Daripa, (2015). If $B_{n}(t)$ is the amplitude of the disturbance with wave number $n$, the growth rate is

[TABLE]

The relationship between equations (27) and (28) is consistent with the comparison of the growth rates in the two different coordinate systems given by equation (24).

We now turn to three-layer flows in which the fluid in the middle layer also has constant viscosity, $\mu_{1}$. This situation has been investigated by Beeson-Jones & Woods, (2015) and Gin & Daripa, (2018), and it has been found that the magnitudes of interfacial disturbances $A_{n}(t)$ and $B_{n}(t)$ are governed by the following system of ODE’s

[TABLE]

where $\mathbf{M}_{1}^{r}(t)$ is the $2\times 2$ matrix with entries given by

[TABLE]

We recall equation (22) which compares the interfacial disturbance of the inner interface in the $r$-coordinates and the $\zeta$-coordinates, and also consider the corresponding equation for the outer interface:

[TABLE]

Equations (29) and (31) gives the matrix equation

[TABLE]

where

[TABLE]

Taking derivatives of (31) and rewriting the resulting equations in matrix form, we obtain

[TABLE]

where

[TABLE]

Combining (32) and (34),

[TABLE]

where $\mathbf{M}_{1}^{\zeta}=\mathbf{R}\mathbf{M}_{1}^{r}\mathbf{R}^{-1}+\mathbf{Q}$. The entries of $\mathbf{M}_{1}^{\zeta}(\tau)$ are

[TABLE]

There are several important things to note from the relationship between the matrices $\mathbf{M}_{1}^{\zeta}$ and $\mathbf{M}_{1}^{r}$:

It was demonstrated by Gin & Daripa, (2018) that $\mathbf{M}_{1}^{r}$ can have complex eigenvalues. However, it is shown in the next section (section 5) that the problem in the $\zeta$-coordinates has real growth rates. This analysis holds even for constant viscosity. Therefore, $\mathbf{M}_{1}^{\zeta}$ has real eigenvalues. 2. 2.

$\mathbf{M}_{1}^{\zeta}$ can also be expressed as $\mathbf{M}_{1}^{\zeta}=\mathbf{R}\left(\mathbf{M}_{1}^{r}+\mathbf{Q}\right)\mathbf{R}^{-1}$. Since a similarity transformation does not change eigenvalues, the eigenvalues of $\mathbf{M}_{1}^{\zeta}$ are the eigenvalues of $\mathbf{M}_{1}^{r}+\mathbf{Q}$ where $\mathbf{Q}$ is a diagonal matrix. Thus it is this diagonal matrix $\mathbf{Q}$ which, when added to $\mathbf{M}_{1}^{r}$, converts the complex eigenvalues to real and leaves real eigenvalues real. 3. 3.

Both $\mathbf{M}_{1}^{r}$ and $\mathbf{M}_{1}^{\zeta}$ have real eigenvalues when $F_{1}$ and $F_{2}$ defined respectively in (16) and (18) have the same sign. Define the matrices:

[TABLE]

then $\mathbf{E}\mathbf{M}_{1}^{r}\mathbf{E}^{-1}$ and $\mathbf{F}\mathbf{M}_{1}^{\zeta}\mathbf{F}^{-1}$ are real symmetric matrices. Therefore, they are similar to real symmetric (i.e. self-adjoint) matrices and have real eigenvalues. However, the argument breaks down when $F_{1}$ and $F_{2}$ have opposite signs because $\mathbf{E}\mathbf{M}_{1}^{r}\mathbf{E}^{-1}$ and $\mathbf{F}\mathbf{M}_{1}^{\zeta}\mathbf{F}^{-1}$ are not symmetric. This shows that the complex eigenvalues of $\mathbf{M}_{1}^{r}$ can only occur when $F_{1}F_{2}\leq 0$.

5 Upper Bounds

To derive an upper bound on the growth rate, we take an inner product of $\eqref{Eigproblem}_{1}$ with $f$. Using integration by parts along with the boundary conditions $\eqref{Eigproblem}_{2}$ and $\eqref{Eigproblem}_{3}$ and solving for $\sigma$ yields

[TABLE]

where

[TABLE]

Note that all terms in (39) are real. Therefore, $\sigma$** is real for all wave numbers.** This is a product of the change of variables from the $r$-coordinates to the $\zeta$-coordinates. It is shown by Gin & Daripa, (2015) that the growth rate can be complex for constant viscosity flows in the $r$-coordinates.

When $\sigma>0$, we may ignore the positive term containing $I_{2}$ in the denominator and get

[TABLE]

We use the following inequality

[TABLE]

which holds for any $N$ if $A_{i}>0$, $B_{i}>0$, and $X_{i}>0$ for all $i=1,...,N$. By using this inequality with $N=3$,

[TABLE]

But

[TABLE]

Therefore,

[TABLE]

Using the definitions of $F_{1}$ and $F_{2}$ given by (16) and (18),

[TABLE]

which is the modal upper bound for a wave with wave number $n$. We can find an absolute upper bound for all wave numbers by taking the maximum of the first two terms over all values of $n$. The absolute upper bound is

[TABLE]

6 Characterization of the Eigenvalues and Eigenfunctions

Using $\lambda=1/\sigma$, the eigenvalue problem (21) can be written as

[TABLE]

Note that $F_{1}$ and $F_{2}$ are positive for small values of $n$ and negative for large values of $n$ (see equations (16) and (18)). From the upper bound (43), we can see that as long as the viscous gradient $\mu^{\prime}(\zeta)$ is not too large, the maximum value of $\sigma$ will occur when $F_{1}$ and $F_{2}$ are positive. For this range of wave numbers, we have the following characterization of the eigenvalues and eigenfunctions.

Theorem 6.1**.**

Let $F_{1}$, $F_{2}$, $Q$, $n$, $\mu_{i}$, $\mu_{o}>0$. Let $\mu(\zeta)$ be a positive, strictly increasing function in $C^{1}([\zeta_{1},1])$. Then the eigenvalue problem (46) has a countably infinite number of real eigenvalues that can be ordered

[TABLE]

with the property that for the corresponding eigenfunctions, $\left\{f_{i}\right\}_{i=0}^{\infty}$, $f_{i}$ has exactly $i$ zeros in the interval $(\zeta_{1},1)$. Additionally, the eigenfunctions are continuous with a continuous derivative.

Proof.

The fact that there are a countably infinite number of real eigenvalues that can be ordered and corresponding eigenfunctions with the prescribed number of zeros is proven by Ince (Ince,, 1956, p. 232-233) in Theorem I and Theorem II using

[TABLE]

The regularity of the eigenfunctions comes from the existence theorem of Ince (Ince,, 1956, p. 73). We saw from equation (39) that $\sigma$ is real for all $n$, and a closer look at each term in (39) shows that if $F_{1},F_{2}>0$ and $\mu(\zeta),\mu^{\prime}(\zeta)>0$, then all terms are positive and $\sigma>0$. ∎

6.1 Self-Adjointness and Expansion Theorem

We rewrite equation (46) as

[TABLE]

This is of the form

[TABLE]

where

[TABLE]

Given the same assumptions as in Theorem 6.1, we have the following theorem from Walter, (1973).

Theorem 6.2**.**

Let $F_{1}$, $F_{2}$, $Q$, $n$, $\mu_{i}$, $\mu_{o}>0$. Let $\mu(\zeta)$ be a positive, strictly increasing function in $C^{1}([\zeta_{1},1])$. Let $p(\zeta)$, $q(\zeta)$, and $r(\zeta)$ be defined by (49). Let

[TABLE]

and define the operator $T$ on $L^{2}_{r}(\zeta_{1},1)$ by

[TABLE]

Define the measure:

[TABLE]

We consider the Hilbert space $H:=L^{2}([\zeta_{1},1];\nu)$. Consider the operator $A$ with domain

[TABLE]

and defined by

[TABLE]

Then $(f,\lambda)$ satisfies (47) if and only if $Af=\lambda f$. $A$ is a self-adjoint operator on $H$ and for any $u\in H$,

[TABLE]

where the $f_{k}$ are the eigenfunctions of $A$.

6.2 Notes on the Assumptions

The above theorem holds assuming

[TABLE]

and $p(\zeta)>0,r(\zeta)>0$ for $\zeta\in[\zeta_{1},1]$. This is satisfied if

[TABLE]

It is also assumed that $F_{1},F_{2}>0$ which holds when

[TABLE]

A different theory will be necessary to consider non-monotonic viscous profiles or large wave numbers.

7 Numerical Results

We now investigate the growth rate of disturbances by numerically computing the eigenvalues of the eigenvalue problem (46). This eigenvalue problem has time-dependent coefficients and boundary conditions which depend on the eigenvalues. Thus the dispersion relation for this problem depends on time. The eigenvalues are computed using a pseudo-spectral Chebyshev method. The eigenvalues $\lambda$ are then inverted to find the growth rates $\sigma$. Recall that for a given wave number $n$, there are infinitely many eigenvalues. In the results that follow, $\sigma$ refers to the maximum over all eigenvalues. $\sigma_{max}$ refers to the maximum over all eigenvalues and over all wave numbers. For consistency, we often use the same parameter values throughout our results. Unless otherwise stated, $\mu_{i}=2$, $\mu_{o}=10$, $T_{1}=T_{2}=1$, and $Q=10$. Therefore, the inner and outer layer fluids have constant viscosity 2 and 10 respectively for all our studies here. The viscous profile of the middle layer fluid, however, is a free variable which can be taken as constant or variable in our studies below. In the rest of the paper, we will characterize the flow by the viscosity of the middle layer.

7.1 Constant vs. Variable Viscosity

We begin by comparing the growth rate of disturbances for a constant viscosity profile with that for a variable viscous profile. In Figure 2(a), the dispersion relations at a fixed time ($\tau=0$) are plotted for four different viscous profiles shown in Figure 2(b). The constant viscosity case is given by the solid (black) line and the viscosity of the middle layer fluid is $\mu=6$. The stability of three-layer constant viscosity flows has been studied extensively by Gin & Daripa, (2015, 2018). Note that there is a maximum growth rate and that short waves are stable due to interfacial tension. For comparison, three linear viscous profiles are considered. The dotted (red) line corresponds to a linear viscous profile with $\mu(R_{1})=5.9$ and $\mu(R_{2})=6.1$, the dashed (blue) line corresponds to $\mu(R_{1})=5$ and $\mu(R_{2})=7$, and the dash-dot (green) line corresponds to $\mu(R_{1})=4$ and $\mu(R_{2})=8$. There are several important features to notice. First, the dispersion relation for each of the variable viscous profiles has a local maximum for a wave number that is similar to the maximum for the case of a constant viscous profile. For profiles with smaller viscous jumps at the interfaces, the local maximum is smaller. Therefore, this local maximum can be attributed to the instability of the interfaces due to the positive viscous jump. The second thing to notice is that short waves are unstable for variable viscous profiles, even when the viscous profile is nearly constant (see the dotted line). As the gradient of the viscous profiles increase, the growth rate of short waves also increases. Therefore, the short wave behavior is dominated by the instability of the middle layer fluid itself due to the viscous gradient. The final observation which can be drawn from Figure 2 is that the maximum growth rate can be smaller for a variable viscous profile than it is for a constant viscous profile with constant viscosity equal to the average of the values of viscosity at the two interfaces in the middle layer of the variable viscosity profile.

7.2 Optimal Profile

In subsection 7.1, it is shown that some particular variable viscous profiles are less unstable (i.e. have a smaller maximum growth rate) than a particular constant viscous profile. This leads to some more general questions: Are there variable viscous profiles that are less unstable than all constant viscous profiles? What is the optimal viscous profile?

The question of the optimal viscous profile is a difficult one so we start by using some simplifying assumptions. First, recall that the viscous profile and the growth rate are both time-dependent. In this section we only consider the growth rate at time $\tau=0$. This is reasonable because it is advantageous to control the instability at early times. Therefore, for the present discussion the term “optimal” refers to the viscous profile that minimizes the maximum growth rate $\sigma_{max}$ at time $\tau=0$. The second simplification is that we first consider only viscous profiles in the middle layer that are linear at time $\tau=0$. Note that linear profiles are uniquely determined by the values $\mu(R_{1})$ and $\mu(R_{2})$. Other types of viscous profiles will be considered later in this section.

Figure 3(a) shows the value of $\sigma_{max}$ for each linear viscous profile such that the viscosity of the middle layer is between $\mu_{i}$ and $\mu_{o}>\mu_{i}$. The optimal viscous profile in this case is $\mu(R_{1})=3.41$ and $\mu(R_{2})=5.09$. Note that all possible constant viscous profiles have been considered as a subset of the set of linear viscous profiles. Therefore, the fact that the optimal profile is not constant shows that variable viscous profiles can be used to reduce the instability of a flow. The dispersion curve for the optimal viscous profile shown in Figure 3 approaches the value indicated by the dotted horizontal line as $n\to\infty$. The two local maxims in this plot have the same value as this limit. This is because the optimal viscous profile is the one which balances the instabilities of the interfaces with the instability of the middle layer.

Next we investigate the optimal viscous profile under several different values of interfacial tension. Plots of $\sigma_{max}$ versus the different linear profiles is given in Figure 4. Figure 4(a) has the smallest value of interfacial tension with $T_{1}=T_{2}=0.25$. The optimal viscous profile has endpoint viscosities of $\mu(R_{1})=3.20$ and $\mu(R_{2})=5.65$. Figure 4(b) uses $T_{1}=T_{2}=1$ and is a repeat of Figure 3(a). As noted above, the optimal profile is $\mu(R_{1})=3.41$ and $\mu(R_{2})=5.09$. Figure 4(c) has the largest values of interfacial tension with $T_{1}=T_{2}=4$. The optimal linear viscous profile is $\mu(R_{1})=3.47$ and $\mu(R_{2})=4.53$. The trend is that larger values of interfacial tension correspond to optimal viscous profiles with a smaller viscous gradient. This is because, as mentioned previously, the optimal viscous profile is the one which balances the instabilities of the interfaces with the instability of the middle layer. A larger value of interfacial tension decreases the instability of the interfaces so the gradient of the middle layer must also decrease in order to match the interfacial instability.

In the case of chemical EOR by polymer flooding, the viscous profile of the middle layer fluid is determined by the concentration of polymer. The use of large quantities of polymer can be expensive so it is useful to explore which viscous profile minimizes the instability (i.e. $\sigma_{\rm max}$) given a fixed amount of total polymer. Assuming that there is a linear relationship between the concentration of polymer and the viscosity, this can be viewed as minimizing the instability for a fixed value of average viscosity. The results of this type of optimization are given in Figure 5. The value of $\sigma_{max}$ for all linear viscous profiles is plotted in Figure 5(a) using the same parameter values as Figure 3(a). For each possible average viscosity between $\mu_{i}$ and $\mu_{o}$, the profile which minimizes $\sigma_{max}$ was found and is marked by an ‘$x$’ in Figure 5(a). The viscous profile can be identified by its slope $a=(\mu(R_{2})-\mu(R_{1}))/(R_{2}-R_{1})$. The slopes of the optimal profiles are plotted versus the average viscosity of the middle layer in Figure 5(b).

When the average viscosity of the middle layer fluid is $\mu_{i}=2$, the optimal viscous profile is constant at $\mu\equiv 2$ (note that this is the only profile considered since we are taking $\mu_{i}\leq\mu(r)\leq\mu_{o}$ for $R_{1}\leq r\leq R_{2}$). Therefore there is no viscous jump at the inner interface and no instability in the layer itself. All of the instability occurs due to the jump in viscosity at the outer interface. As the average viscosity increases from there, the jump at the inner interface of the optimal viscous profile remains zero while the viscous gradient increases in order to decrease the viscous jump at the outer interface. Eventually, as the average viscosity nears $\mu=3$, a point is reached in which the gradient stops increasing as illustrated by the flat portion of the graph in Figure 5(b). During this time, the viscous jump at the inner interface increases while the viscous jump at the outer interface decreases. The point at which the slope begins to increase again corresponds to the optimal viscous profile over all values of average viscosity. After this point, the addition of polymer to increase viscosity would be detrimental to the stability of the system.

To this point, only linear viscous profiles have been considered, but many other viscous profiles can also be used. The optimization procedure used for linear viscous profiles in Figure 3(a) has been repeated for several other types of viscous profiles in Figure 6. In addition to a viscous profile which is linear at time $\tau=0$, viscous profiles which are initially exponential and logarithmic are also considered. Also, recall that the viscous profile changes with time in the $r$-coordinates. Therefore, it may be useful to consider viscous profiles in the $\zeta$-coordinate system because they will remain fixed in time. Therefore, we also consider viscous profiles which are linear and exponential with respect to $\zeta$. A profile which is linear with respect to $\zeta$ is quadratic with respect to $r$ and a profile which is exponential in $\zeta$ is proportional to $e^{r^{2}}$. Figure 6(a) shows the dispersion relations of the optimal viscous profiles of each type. The profile which is exponential in $r$ is the most unstable while the profile which is exponential in $\zeta$ is the least unstable. The corresponding optimal viscous profiles are plotted in Figure 6(b). Notice that the profile which is least unstable, the one which is exponential in $\zeta$, has both the smallest value of $\mu(R_{1})$ and the largest value of $\mu(R_{2})$. Therefore, it has the smallest viscous jumps at the interfaces.

7.3 Time dependence

In the previous sections, we considered the growth rate only at time $\tau=0$. However, it is also important to understand how the growth rate changes with time. As time increases and the interfaces move outward, there are several physical factors at play. The curvature of the interfaces decreases which works to destabilize the flow while the velocity of the interfaces decreases which works to stabilize the flow. Gin & Daripa, (2015) studied constant viscosity flows and found that the above two competing effects lead to non-monotonic behavior of the maximum growth rate. However, the growth rate of constant viscosity flows behaves differently in the $\zeta$-coordinate system. We illustrate this first with some analytical results.

For two-layer constant viscosity flows, the growth rate in the $\zeta$-coordinates (see equation (27)) is

[TABLE]

where $R(\tau)$ is the radius of the interface. By taking a derivative with respect to $n$ and setting equal to zero, the most dangerous wave number is

$n_{max}=\sqrt{QR/(6\pi T)(\mu_{o}-\mu_{i})+1/3}$. Plugging this into (27) gives

[TABLE]

Taking the derivative with respect to $R$ gives

[TABLE]

If $\mu_{o}>\mu_{i}$ then this expression is negative for all $R$. Therefore, $\sigma_{max}$ is a strictly decreasing function of $R$ and hence time $\tau$ in the $\zeta$-coordinates.

For three-layer variable viscosity flows, there is an additional factor which affects the stability. The interfaces get closer together which makes the variable viscous profile steeper and works to destabilize the flow. Despite this fact, the numerical results that follow show that $\sigma_{max}$ is a decreasing function of time. This can be illustrated by the upper bound given by (45). The first two terms are strictly decreasing functions of $R_{1}$ and $R_{2}$ (and therefore of $\tau$) while the third term is independent of time. Therefore, the upper bound is a decreasing function of time. In Figure 7, the dispersion relation is plotted at several different times for a typical variable viscosity flow. Initially, the inner interface is at $R_{1}=20$. As time increases, it moves outward. Note that the maximum value of $\sigma$ decreases with time. However, the difference is more pronounced near the maximum value, which is mostly affected by the stability of the interfaces, than for short waves which are mostly affected by the layer instability.

In order to shed more light on this time-dependent behavior, we investigate how $\sigma_{max}$ and the most dangerous wave number $n_{max}$ evolve in time. Figure 8(a) is a plot of $\sigma_{max}$ versus $R_{1}$, and Figure 8(b) is a plot of $n_{max}$ versus $R_{1}$. Notice first that $\sigma_{max}$ is a decreasing function of time and $n_{max}$ is an increasing function of time. The fact that $n_{max}$ increases with time is a well-known fact for constant viscosity flows (Cardoso & Woods,, 1995). Also observe that there is a critical value $R_{1}^{*}$ such that for $R_{1}>R_{1}^{*}$, $\sigma_{max}$ is constant and $n_{max}$ is infinite. This is the point at which the layer instability comes to dominate the flow. For $R_{1}<R_{1}^{*}$, the instability of the interfaces dominates and the behavior or $\sigma_{max}$ and $n_{max}$ is similar to what happens for constant viscosity flow. For $R_{1}>R_{1}^{*}$, the layer is more unstable than the interfaces, and therefore the short wave instability dominates and $\sigma_{max}=\displaystyle{\lim_{n\to\infty}\sigma(n)}$ . This behavior mimics what we see from the upper bound (45) in which the two terms related to the interfaces are decreasing functions while the term related to the layer instability is constant.

7.4 Validation of QSSA

In Section 2, we invoke the QSSA which assumes that the basic solution changes slowly in comparison to the growth of the disturbances. We now examine the validity of that assumption. Consider a two-layer constant viscosity flow in which the growth rate of a disturbance with wave number $n$ is given by equation (27). The maximum growth rate over all wave numbers can be written as

[TABLE]

The expression for the position of the interface is $R(\tau)=\sqrt{Q\tau/\pi+R(0)^{2}}$. In particular, $R\propto\sqrt{\tau}$. Therefore, $\sigma_{max}\propto\tau^{-3/4}$ for $\tau\gg 1$. In comparison, the interfacial position of the basic solution changes at a rate

[TABLE]

Therefore, for large $\tau$ the disturbances will grow faster than the basic solution. For three-layer variable viscosity flow, the QSSA has a more solid foundation. Recall from the previous subsection that the interfacial instability dominates at early times, but the instability of the middle layer dominates at later times. The upper bound (45) demonstrates that the layer instability is bounded by a constant term that depends only on $\mu^{\prime}(\zeta)$. This is further validated by the region of Figure 8(a) in which $\sigma_{max}$ is constant. Therefore, for variable viscosity flows the interfaces will be moving very slowly at later times while the growth of disturbances remains constant.

Figure 9 is a numerical comparison of the growth rate of disturbances $\sigma_{max}$ and the rates of change of each individual interface of the base flow. The parameters used are the same as in Figure 8 and therefore the solid curve is the same as the curve in Figure 8(a). Notice that the growth rate of the disturbances is always greater than the rate of change of the interfaces, but that this is especially true for later times. The late time behavior will be true even if the interfaces are stabilized by very large interfacial tension.

7.5 Variable Injection Rate

Recent studies by Beeson-Jones & Woods, (2015) and Gin & Daripa, (2018) explore the idea of using a variable injection rate $Q(t)$ to stabilize multi-layer constant viscosity flows. In these works, the maximum injection rate which results in a stable flow is investigated. Unfortunately, there is no injection rate which stabilizes a variable viscosity flow because short waves are always unstable. However, as an analogy, we can find the maximum injection rate that keeps the growth rate under a certain threshold. Figure 10 shows the maximum injection rate such that the maximum growth rate is below 0.001 for a certain constant viscosity flow and a certain variable viscosity flow. The constant viscosity flow has a viscosity of $\mu=6$ in the middle layer while the variable viscosity flow has a linear viscous profile with $\mu(R_{1})=5$ and $\mu(R_{2})=7$. The variable viscosity flow allows for the fluid to be injected more quickly while maintaining the same level of instability.

8 Conclusions

The stability of three-layer radial porous media flows with variable viscosity is an important issue in many applications. This work is the first to address this topic. First, the linear stability problem is formulated resulting in an eigenvalue problem with time-dependent coefficients and eigenvalue-dependent boundary conditions. This derivation requires an appropriate change of variables that fixes the positions of the interfaces and the viscous profile of the middle layer fluid. Several important analytical aspects of the eigenvalue problem are studied. First, upper bounds on the spectrum are derived using a variational approach. Then it is shown that for a certain bandwidth of wave numbers, there is a countably infinite set of positive eigenvalues, and the corresponding eigenfunctions are complete in a certain $L^{2}$ space.

The eigenvalues have been computed numerically in order to investigate the effect of various parameters on the stability of variable viscosity flows. The following are some of the key numerical results: (i) Variable viscosity flows can reduce the maximum growth rate by reducing the jumps in viscosity at the interfaces, but short waves become unstable; (ii) Short wave instability is dominated by the viscous gradient in the layer while long and intermediate wavelengths are dominated by the instability of the interfaces; (iii) The optimal viscous profile is the one which balances the interfacial instability with the instability of the layer; (iv) increasing interfacial tension decreases the viscous gradient of the optimal viscous profile; (v) A viscous profile which is exponential with respect to $\zeta$ is optimal among the types of profiles considered; (vi) $\sigma_{max}$ is a decreasing function of time. This is due to the instability of the interfaces decreasing with time while the layer instability remains relatively unchanged; and (vii) Variable viscosity flows allow for faster injection without making the flow more unstable.

Acknowledgment

This work has been supported in part by the U.S. National Science Foundation grant DMS-1522782.