1. Introduction
The parallel co-flow of two fluids occurs in many industrial, biological and environmental processes. It is often important to understand the interfacial instability and develop strategies to control it. Frequently, these flows occur in thin channels whose thickness varies in the cross-flow direction. Examples include: the flow of cement and drilling fluid within a casing pipe of a subsurface well, where the intermingling of cement and mud can produce poorly sealed wells with the attendant risks of leakage; the flow of coatings in the corner region along the line of intersection between two planes, where the displacement of air limits the formation of non-coated zones (Weislogel & Lichter Reference Weislogel and Lichter1996); flows of reactants in microfluidic channels (Sauer Reference Sauer1987; Huang et al. Reference Huang, Yang, Su, Li, Hu, Li, Pan, Ren, Li and Song2018); and the displacement of water by 
$\textrm {CO}_{2}$ in permeable channels used for 
$\textrm {CO}_{2}$ sequestration, where intermingling may enhance the efficiency of the sequestration (Woods & Mingotti Reference Woods and Mingotti2016). For the last example of 
$\textrm {CO}_{2}$ sequestration, the pore-scale dynamics, which is controlled by capillary effects and the interpore geometry, is not fully understood but can have a significant effect on macroscale mechanisms such as the flux and the residual trapping of the 
$\textrm {CO}_{2}$, which are key to estimating storage efficiency (Zhao et al. Reference Zhao2019; Benham, Bickle & Neufeld Reference Benham, Bickle and Neufeld2021). If the initial flow involves the displacement of one fluid by a second along the channel, the displacing fluid will migrate along the wide part of the channel, stretching out the interface; at long times, the interface is approximately directed along the channel, and the flow evolves to the co-flowing geometry of the present problem (cf. Woods & Mingotti Reference Woods and Mingotti2016; Mortimer & Woods Reference Mortimer and Woods2021).
We investigate how the stability of the interface between the two fluids in a thin cell with vertically varying gap width is controlled by cross-layer buoyancy and capillary effects (figure 1). We assume that inertia plays a negligible role in the base flow. It has been shown that when inertia plays a significant role and the two fluids in the Hele-Shaw cell have significantly different viscosities, the shear-controlled Kelvin–Helmholtz instability may occur as has been observed experimentally (Zeybek & Yortsos Reference Zeybek and Yortsos1992; Gondret & Rabaud Reference Gondret and Rabaud1997; Rabaud & Moisy Reference Rabaud and Moisy2020). In wider channels, it has been shown that even at zero Reynolds number the vertical shear associated with the no-slip boundaries at the top and bottom can give rise to interfacial instabilities in the co-flow of two fluids of different viscosity (Yih Reference Yih1967). Similar behaviour can occur in two-layer gravity-driven flow (Loewenherz & Lawrence Reference Loewenherz and Lawrence1989).
(a) Schematic of the set-up. Large arrows indicate the flow direction. (b) Cross-section perpendicular to the 
$x$ direction. (c) Cross-section in the 
$x$ direction.
In the present work, we consider a laterally thin cell in which the vertically varying velocity arises from variations in the cell width. The stability is primarily controlled by the density difference between the fluids and surface tension. It is well established that the along-flow component of surface tension stabilises larger wavenumbers. The combination of surface tension and the cross-cell variation in thickness introduces a new (de)stabilising process for small wavenumbers in the case that the (non-)wetting fluid is in the thinner part of the channel. This effect may complement or suppress the effect of a density difference between the two fluids on the stability of the interface.
The impact of variations in the surface tension associated with variations in the channel width have been explored in detail for the related problem in which an input fluid displaces an ambient fluid in a cell whose width varies in the direction of flow (Homsy Reference Homsy1987; Al-Housseiny, Tsai & Stone Reference Al-Housseiny, Tsai and Stone2012; Dias & Miranda Reference Dias and Miranda2013; Grenfell-Shaw & Woods Reference Grenfell-Shaw and Woods2017). These studies have identified that the effect of cross-cell curvature can complement or suppress the classical Saffman–Taylor instability (Saffman & Taylor Reference Saffman and Taylor1958).
2. Formulation
The flow and the cell geometry are illustrated in figure 1. The cell occupies 
$0< y< H$ and has width 
$b(y)$ which varies in the vertical direction,
\begin{equation} b(y) = b_0 + \alpha y, \end{equation}
where 
$\alpha$ represents the inclination of the cell walls, which may be positive or negative, and 
$\alpha >-b_0/H$, so that the cell width is non-negative. Flow is driven in both fluids in the 
$x$ direction by a background pressure gradient with magnitude 
$G$. For relatively slow flows, we can apply the lubrication approximation, which is to say that the leading-order velocities are independent of 
$x$. Under this assumption, the momentum and continuity equations for the gap-averaged velocity in each fluid take the form (equation (11) of Gondret & Rabaud Reference Gondret and Rabaud1997)
\begin{equation} \frac{\partial \bar{\boldsymbol{u}}}{\partial t} + \bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{u}}={-}\frac{1}{\rho} \boldsymbol{\nabla} p -\frac{12 \mu}{\rho b^2} \bar{\boldsymbol{u}} - g \boldsymbol{e}_y, \quad \boldsymbol{\nabla} \boldsymbol{\cdot}(b \bar{\boldsymbol{u}})=0, \end{equation}
where 
$\boldsymbol {\nabla }=(\partial /\partial x,\partial /\partial y)$, 
$\boldsymbol {e}_y$ is the unit vector in the 
$y$ direction and 
$\mu$ and 
$\rho$ are the fluid viscosity and density, respectively. Also, 
$p$ is the pressure, 
$g$ represents gravity, which acts in the negative 
$y$ direction, and 
$\bar {\boldsymbol {u}}=(\bar {u},\bar {v})$ is the width-averaged velocity in the 
$x$ and 
$y$ directions. The boundary conditions are no-flux at the top and bottom of the channel, 
$\bar {v}=0$, whilst at the fluid–fluid interface, 
$y=y_I$, the velocity in each fluid satisfies the kinematic boundary condition
\begin{equation} \frac{\partial y_I}{\partial t} + \bar{u} \frac{\partial y_I}{\partial x} = \bar{v}. \end{equation}
In addition, there is a pressure jump at the interface associated with its curvature, 
$\kappa =\nabla ^2 y_I$, given by
\begin{equation} \Delta p = \gamma \kappa, \end{equation}
where 
$\gamma$ represents surface tension. The unperturbed steady base flow is given by
\begin{equation} \bar{u} = U_0(y)= \frac{b(y)^2 G}{12 \mu}, \quad \bar{v}=0, \quad p=P_0(x,y)={-}\rho g y -G x + \text{const}. \end{equation}
The fluids are immiscible and the location of the fluid–fluid interface, 
$y_I=h$, is a constant for the case of steady flow, which depends on the channel angle, the relative flux in each layer and the viscosity ratio. Although the curvature of the interface in the along-flow direction vanishes since 
$y_I$ is independent of 
$x$, there is curvature in the cross-channel direction owing to the contact angle at the wall and the varying channel width (figure 1b). Hence, there is a pressure jump at the interface, which is independent of 
$x$, and the constants in 
$P_0$ in the base flow are different in the two fluids.
We consider perturbations to the interface and steady base flow of the form
$$\begin{gather} y_I= h+ \zeta \exp({\textrm{i}(kx-\omega t)}), \end{gather}$$
$$\begin{gather}(\bar{u},\bar{v}) =(U_0(y),0) +(u(y),v(y)) \exp({\textrm{i}(kx-\omega t)}), \end{gather}$$
$$\begin{gather}p = P_0(x,y) + p(y) \exp({\textrm{i}(kx-\omega t)}), \end{gather}$$
where 
$\zeta$, 
$u(y)$, 
$v(y)$ and 
$p(y)$ are assumed to be small. We seek to determine the stability of such perturbations. The linearised governing equations in each fluid are
$$\begin{gather} - \textrm{i} \omega u + \textrm{i} k U_0 u + v U_0' = \frac{-\textrm{i}kp}{\rho} - \frac{12 \mu u}{\rho b^2}, \end{gather}$$
$$\begin{gather}- \textrm{i} \omega v + \textrm{i} k U_0 v ={-}\frac{p'}{\rho} - \frac{12 \mu v}{\rho b^2}, \end{gather}$$
$$\begin{gather}\textrm{i}kub+ (vb)' =0, \end{gather}$$
where a prime (
$'$) denotes differentiation with respect to 
$y$. We eliminate 
$p'$ to obtain
\begin{equation} -\textrm{i} \omega u' + \textrm{i}k (U_0 u)' +(v U_0')' - k \omega v + k^2 U_0 v= \frac{12 \textrm{i} k \mu v}{\rho b^2} - \frac{12 \mu (u/b^2)'}{\rho}. \end{equation}Also, continuity yields
\begin{equation} u = \frac{\textrm{i} (vb)'}{k b}=\frac{\textrm{i}}{k} \left(v'+\frac{v \alpha}{b} \right), \end{equation}
which can be used to eliminate 
$u(y)$ from (2.12) and obtain a differential equation for 
$v(y)$,
\begin{equation} A(y) v'' + B(y) v' + C(y) v = 0, \end{equation}where the coefficients are defined below for the dimensionless analogue.
2.1. Non-dimensionalisation
To scale the system, we use the tank height, 
$H$, as the length-scale, and the time-scale is 
$T=\mu _2/(GH)$. The pressure scale is 
$GH$. We write
\begin{equation} (\hat{x}, \hat{y}) = (x, y)/H, \quad \hat{k}= H k, \quad \hat{w} = T w, \quad \hat{b}(\hat{y}) = b(y)/H= \widehat{b_0} + \alpha \hat{y}, \end{equation}
with 
$\hat {b}_0=b_0/H$, 
$\hat {h}=h/H$ and 
$\hat {\zeta }=\zeta /H$. The upper fluid is labelled fluid 2, whilst the lower is labelled fluid 1 (figure 1). Henceforth, all quantities are dimensionless unless stated otherwise, and we discard the hat notation. The dimensionless equation in fluid 
$j=1,2$ becomes
\begin{equation} A^{(\,j)}(y) v^{(\,j)}_{yy}+ B^{(\,j)}(y) v^{(\,j)}_{y} + C^{(\,j)}(y) v^{(\,j)} = 0, \end{equation}with coefficients
$$\begin{gather} A^{(\,j)}= \frac{\omega b(y)^4}{k} - \frac{M_j b(y)^6}{12} + \frac{\textrm{i} b(y)^2 R_j\mathcal{A}}{k}, \end{gather}$$
$$\begin{gather}B^{(\,j)} = \frac{\omega \alpha b(y)^3}{k} - \frac{M_j\alpha b(y)^5}{12} - \frac{ \textrm{i} \alpha b(y) R_j\mathcal{A}}{k}, \end{gather}$$
$$\begin{gather}C^{(\,j)} = \frac{-\omega \alpha^2 b(y)^2}{k} +\frac{M_j \alpha^2 b(y)^4}{12} - k \omega b(y)^4 + \frac{M_j k^2 b(y)^6}{12} - \textrm{i} k b(y)^2 R_j\mathcal{A} - \frac{3 \textrm{i} \alpha^2 R_j\mathcal{A}}{k}, \end{gather}$$
where 
$M_1=M$, 
$M_2=1$ and 
$R_1=R$, 
$R_2=M$, and we have introduced the dimensionless parameters
\begin{equation} M=\frac{\mu_2}{\mu_1},\quad R = \frac{\rho_2}{\rho_1}, \quad \mathcal{A}=\frac{12 \mu_1 \mu_2}{\rho_2 H^3 G}, \end{equation}
which respectively represent the viscosity ratio, the density ratio and the importance of viscous drag relative to inertia; 
$\mathcal {A}$ is inversely proportional to a Reynolds number.
2.2. Boundary conditions
The perturbed no-flux boundary conditions are
\begin{equation} v^{(1)}(0) = v^{(2)}(1) = 0. \end{equation}The kinematic boundary conditions (2.3) at the fluid–fluid interface become
\begin{equation} v^{(1)}(h) = \textrm{i} \zeta \left[\frac{b(h)^2 M}{12} k- \omega \right], \quad v^{(2)}(h) = \textrm{i} \zeta \left[\frac{b(h)^2 }{12} k- \omega \right]. \end{equation}
The dynamic boundary condition at the interface accounts for a jump in pressure associated with surface tension and curvature. This comprises two contributions: the along-channel curvature and the cross-channel curvature. The former is proportional to the second derivative of the interface 
$y_I$ in the 
$x$ direction, which furnishes a term proportional to 
$k^2$. Treating the cross-channel curvature requires more care. The contact angle, 
$\tilde {\theta }$, is defined as the angle between fluid 1 and the channel wall (figure 1b). In a cell with inclined walls, the radius of curvature is adjusted from the case of a parallel-sided cell and we define an effective contact angle, 
$\theta = \tilde {\theta }-\phi$, where 
$\tan \phi = \alpha /2$ (Park & Homsy Reference Park and Homsy1984; Romero & Yost Reference Romero and Yost1996; Grenfell-Shaw & Woods Reference Grenfell-Shaw and Woods2017). The discontinuity in the perturbed pressures at the interface is then given by
\begin{equation} p_2(h) - p_1(h) = C \zeta \left[ R-1 - {\textit{Bo}}^{{-}1}\left(k^2 + \frac{2 \alpha \cos \theta}{b(h)^2} \right) \right], \end{equation}where we have introduced the following dimensionless groups:
\begin{equation} C=\frac{\rho_1 g}{G}, \quad {\textit{Bo}} = \frac{g\rho_1 H^2}{\gamma}. \end{equation}
We use the dimensionless analogues of (2.9) and (2.13) to obtain the dimensionless pressures in terms of the vertical velocity, 
$v$, in each fluid,
\begin{equation} p_j = \frac{12 \textrm{i} \omega}{k^2 R_j M_j \mathcal{A}}\frac{(v^{(\,j)} b)'}{b}-\frac{\textrm{i}}{k}\frac{b (v^{(\,j)} b)'}{\mathcal{A} R_j} + \frac{2\textrm{i} \alpha b v^{(\,j)}}{k\mathcal{A}R_j}- \frac{12}{M_j k^2} \frac{(v^{(\,j)}b)'}{b^3}. \end{equation}We substitute the pressures into the dynamic boundary condition (2.23) to obtain
\begin{align} &\alpha v^{(1)} \left(\frac{\omega}{kb}+ \frac{Mb}{12}+\frac{\textrm{i} R\mathcal{A}}{kb^3}\right) + \frac{\mathrm{d} v^{(1)}}{\mathrm{d}\kern0.05em y} \left(\frac{\omega}{k} - \frac{M b^2}{12}+\frac{\textrm{i} R\mathcal{A}}{kb^2}\right) \nonumber\\ &\qquad -R \alpha v^{(2)} \left(\frac{\omega}{kb}+ \frac{b}{12}+\frac{\textrm{i} M\mathcal{A}}{kb^3}\right) - R \frac{\mathrm{d} v^{(2)}}{\mathrm{d} y} \left(\frac{\omega}{k} - \frac{b^2}{12}+\frac{\textrm{i} M\mathcal{A}}{kb^2}\right)\nonumber\\ &\quad=\frac{k R M\mathcal{A}\textrm{i}}{12} C \zeta \left[ R-1 - {\textit{Bo}}^{{-}1}\left( k^2 + \frac{2 \alpha \cos \theta}{b^2} \right) \right] , \end{align}
where the suppressed argument of 
$v^{(1)}$, 
$v^{(2)}$, 
$\mathrm {d} v^{(1)}/\mathrm {d} y$, 
$\mathrm {d} v^{(2)}/\mathrm {d} y$ and 
$b$ is 
$y=h$.
3. Solution method
The system for 
$v(y)$ comprises two second-order linear ordinary differential equations in 
$0< y< h$ and 
$h< y< H$ (2.16) and four boundary conditions: no-flux at the top and bottom boundaries (2.21) and the kinematic and dynamic boundary conditions (2.22a,b), (2.26) at the interface. We note that the two equations for the kinematic condition (2.22a,b) can be used to eliminate 
$\zeta$ from the problem. To solve this system, we first simplify the problem by writing 
$\bar {b}(y)=kb(y)/\alpha$ and obtain the following equation for 
$v^{(\,j)}(\bar {b})$:
\begin{equation} A^{(\,j)}(y) v^{(\,j)}_{\bar{b}\bar{b}}+ B^{(\,j)}(y) v^{(\,j)}_{\bar{b}} + C^{(\,j)}(y) v^{(\,j)} = 0, \end{equation}with
$$\begin{gather} A^{(\,j)}= a_6^{(\,j)} \bar{b}^6 + a_4^{(\,j)} \bar{b}^4 + \bar{b}^2, \end{gather}$$
$$\begin{gather}B^{(\,j)} = a_6^{(\,j)}\bar{b}^5 + a_4^{(\,j)} \bar{b}^3 -\bar{b}, \end{gather}$$
$$\begin{gather}C^{(\,j)} ={-}a_6^{(\,j)}(\bar{b}^6 +\bar{b}^4)-a_4^{(\,j)} (\bar{b}^4 + \bar{b}^2) - \bar{b}^2 -3, \end{gather}$$where for each fluid we have introduced
\begin{equation} a_6^{(\,j)}= \frac{-M_j \alpha^4}{12\textrm{i}R_j \mathcal{A}k^3}, \quad a_4^{(\,j)}=\frac{\omega \alpha^2}{\textrm{i}R_j \mathcal{A}k^2}. \end{equation} The general solution for 
$v$ in each fluid is given by a linear combination of two independent power series, 
$\varPhi ^{(\,j)}(\bar {b})$ and 
$\varPsi ^{(\,j)}(\bar {b})$, whose coefficients are determined by Frobenius’ method (given in Appendix A). The velocities are (
$j=1,2$)
\begin{equation} v^{(\,j)}(y) =
c^{(\,j)} \left(
\frac{\varPhi^{(\,j)}(\bar{b}(y))}{\varPhi^{(\,j)}(\bar{b}(\,j-1))}
-
\frac{\varPsi^{(\,j)}(\bar{b}(y))}{\varPsi^{(\,j)}(\bar{b}(\,j-1))}\right),
\end{equation}
where 
$c^{(\,j)}$ are constants and we have used the no-flux boundary conditions at the base and the top (
$y=j-1$). The dynamic boundary condition may be written as follows (at 
$y=h$):
\begin{align} &\frac{(v^{(1)} \bar{b})_{\bar{b}}}{\bar{b}^3} + a_4^{(1)} \frac{(v^{(1)} \bar{b})_{\bar{b}}}{\bar{b}} + a_6^{(1)}\left( \frac{v^{(1)}}{\bar{b}}\right)_{\bar{b}} \bar{b}^3 \end{align}
\begin{align} &\qquad -M \left[ \frac{(v^{(2)} \bar{b})_{\bar{b}}}{\bar{b}^3} + a_4^{(2)} \frac{(v^{(2)} \bar{b})_{\bar{b}}}{\bar{b}} + a_6^{(2)}\left( \frac{v^{(2)}}{\bar{b}}\right)_{\bar{b}} \bar{b}^3 \right] \end{align}
\begin{align} &\quad=\frac{\textrm{i}M \alpha^2 v^{(1)}}{12k \omega} C \left[ R-1 - {\textit{Bo}}^{{-}1}\left( k^2 + \frac{2 \alpha \cos \theta}{b^2} \right) \right] \left(\frac{a_6^{(1)}}{a_4^{(1)}} \bar{b}^2 +1 \right)^{{-}1}. \end{align}The kinematic boundary condition may be written as
\begin{equation} v^{(1)} \left(\frac{a_6^{(2)}}{a_4^{(2)}} \bar{b}^2 +1 \right)=v^{(2)} \left(\frac{a_6^{(1)}}{a_4^{(1)}} \bar{b}^2 +1 \right). \end{equation}Combining the velocities with these two boundary conditions furnishes the following dispersion relation:
\begin{equation} D^{(1)}E^{(1)}-M E^{(2)} D^{(2)} -S = 0, \end{equation}where
$$\begin{gather} E^{(\,j)}=\frac{t_j}{\bar{b}^2} +\frac{1}{\bar{b}^3} +a_4^{(\,j)} \left(t_j+\frac{1}{\bar{b}} \right) +a_6^{(\,j)} \left(t_j \bar{b}^2-\bar{b} \right), \end{gather}$$
$$\begin{gather}D^{(\,j)} = \frac{a_6^{(\,j)}}{a_4^{(\,j)}} \bar{b}^2 +1, \end{gather}$$
$$\begin{gather}S = \frac{\textrm{i}M\alpha^2}{12k \omega}C \left[R-1-{\textit{Bo}}^{{-}1}(k^2+2\alpha\cos \theta/b(h)^2) \right], \end{gather}$$
$$\begin{gather}t_1 = \frac{\varPhi^{(1)}_{\bar{b}}(\bar{b}(h))/\varPhi^{(1)}(\bar{b}(0))-\varPsi^{(1)}_{\bar{b}}(\bar{b}(h))/ \varPsi^{(1)}(\bar{b}(0))}{\varPhi^{(1)}(\bar{b}(h))/\varPhi^{(1)}(\bar{b}(0))-\varPsi^{(1)}(\bar{b}(h))/ \varPsi^{(1)}(\bar{b}(0))}, \end{gather}$$
$$\begin{gather}t_2 = \frac{\varPhi^{(2)}_{\bar{b}}(\bar{b}(h))/\varPhi^{(2)}(\bar{b}(1))-\varPsi^{(2)}_{\bar{b}} (\bar{b}(h))/\varPsi^{(2)}(\bar{b}(1))}{\varPhi^{(2)}(\bar{b}(h))/\varPhi^{(2)}(\bar{b}(1))- \varPsi^{(2)}(\bar{b}(h))/\varPsi^{(2)}(\bar{b}(1))}, \end{gather}$$
with 
$\textrm {i}$ denoting the imaginary unit. For any value of 
$k$ and the dimensionless parameters, we may obtain a solution for the growth rate, 
$\omega$, that satisfies (3.11) (e.g. figure 2a).
(a) Growth rate, 
$\omega _I$, as a function of wavenumber 
$k$ in the case of equal-density fluids (
$R=1$). The curves are calculated using the method in § 3. The critical wavenumbers predicted by (4.3) are shown as crosses. We use 
${\textit {Bo}}=5$, 
$h=0.5$, 
$b_0=0.3$, 
$M=2$, 
$\mathcal {A}=1$ and 
$C=1$. (b) Schematic corresponding to the blue curve.
4. Analysis
The terms in the square brackets in 
$S$ (3.14) correspond to the pressure jump at the interface, and we define
\begin{equation} J=R-1-{\textit{Bo}}^{{-}1}[k^2+2\alpha\cos \theta/b(h)^2]. \end{equation}
In general, 
$J>0$ is associated with instability and 
$J<0$ is associated with stability. The first term, 
$R-1$, represents the density difference between the fluids. It stabilises the interface for 
$R<1$ and destabilises it for 
$R>1$. The term 
$-{\textit {Bo}}^{-1}k^2$ stabilises the interface; it arises from surface tension suppressing the curvature in the along-channel (
$x$) direction. The final term, 
$-2{\textit {Bo}}^{-1}\alpha \cos \theta /b(h)^2$, is associated with surface tension acting on the curvature across the thickness of the cell. It drives or suppresses an instability depending on whether the wetting or non-wetting fluid is in the thinner part of the channel. This corresponds to the sign of 
$\alpha$ and the sign of 
$\cos \theta$. To interpret the instability, we consider the simpler cases of equal density in § 4.1 and parallel walls in § 4.2 before returning to the full problem in § 4.3.
4.1. Equal density (
$R=1$)
In the case of equal-density fluids, 
$R=1$, the pressure jump reduces to
\begin{equation} J={-}{\textit{Bo}}^{{-}1}[ k^2 + 2\alpha\cos \theta/b(h)^2 ]. \end{equation}
For small wavenumber, 
$k$, the cross-cell surface tension term controls the stability as demonstrated in figure 2(a). The red and yellow curves correspond to the non-wetting fluid occupying the thicker part of the channel, and the system is stable, as this fluid will remain in the thicker part of the channel. The blue and purple curves represent the converse situation, in which the non-wetting fluid is in the thinner part of the channel and will move to the wider side of the channel, leading to an instability (see figure 2b). For larger wavenumbers, the along-channel term stabilises the interface. The along-channel and cross-channel curvature terms balance (
${\textit {Bo}}^{-1}k^2 + 2{\textit {Bo}}^{-1}\alpha \cos \theta /b(h)^2=0$) at a critical wavenumber,
\begin{equation} k_c = \frac{\sqrt{-2 \alpha \cos \theta}}{b(h)}. \end{equation}
For 
$k>k_c$, we anticipate that the system is stable. The critical wavenumbers are shown by crosses for the two unstable set-ups in figure 2, demonstrating good agreement with the predictions from § 3. The small discrepancy between the prediction of 
$J$ (4.3) and the numerical results arises because of the physical effects, such as inertia, that are incorporated in the numerics but are neglected when using 
$J$ as an approximation of the stability criterion.
In many settings, it is important to understand how the instability depends on the flux in each layer. To analyse this, we calculate the relative flux, 
$\mathcal {Q}$, of the top to the bottom layer,
\begin{equation} \mathcal{Q} = \frac{Q_2}{Q_1} = M^{{-}1} \frac{ (b_0 +\alpha)^3 -(b_0 +\alpha h)^3}{(b_0 +\alpha h)^3-b_0^3}. \end{equation}
The quantity 
$M \mathcal {Q} = \mu _2 Q_2 / (\mu _1 Q_1)$ depends only on 
$\alpha$, 
$h$ and 
$b_0$. We consider channels of fixed dimensionless area, so that
\begin{equation} \int_0^1 b(y) \, \mathrm{d} y = b_0 + \alpha/2 = \text{const}. \end{equation}
For a given channel area and relative flux 
$\mathcal {Q}$, we can calculate the interface height 
$h$ as a function of the channel angle 
$\alpha$ (see figure 3(a) for the case with dimensionless area 
$0.2$). We can also calculate the critical wavenumber, 
$k_c$, by using (4.3) (see figure 3(b) for the case 
$\theta =3{\rm \pi} /4$).
(a) The interface height, 
$h$, as a function of the relative flux and viscosity and the channel angle, 
$\alpha$, for a fixed dimensionless channel area of 
$0.2$. (b) The corresponding critical wavenumber, 
$k_c$, above which the system is stable according to (4.3) for 
$\theta =3{\rm \pi} /4$. Note that the system is stable for all 
$k$ for 
$\alpha \leqslant 0$.
4.2. Parallel cell walls ($\alpha =0$)
In the case that the cell walls are parallel,
\begin{equation} J=R-1-{\textit{Bo}}^{{-}1}k^2. \end{equation}
For small 
$k$, the stability is controlled solely by the density ratio, 
$R$. In the case that 
$R<1$, the system is stabilised by the density difference and there are no destabilising effects (Gondret & Rabaud Reference Gondret and Rabaud1997). For 
$R>1$, the Rayleigh–Taylor instability is stabilised for large 
$k$ owing to the along-cell surface tension. Neutral stability is given by
\begin{equation} k_c = \sqrt{{\textit{Bo}}(R-1)}. \end{equation}
Figure 4 shows the growth rate as a function of wavenumber (obtained in § 3). For 
$R>1$, the critical wavenumber prediction is indicated by crosses, showing good agreement.
Growth rate, 
$\omega _I$, as a function of wavenumber, 
$k$, in the case of parallel cell walls (
$\alpha =0$) calculated using the method of § 3. The density ratio is 
$R=0.4$, 
$0.8$, 
$1.6$, 
$3.2$, 
$6.4$. The crosses correspond to the critical wavenumber for neutral stability for 
$R>1$ given by (4.7). We use 
$h=0.5$, 
$b_0=0.3$, 
$M=2$, 
$\mathcal {A}=1$, 
$C=1$ and 
${\textit {Bo}}=1$.
4.3. Competition between density difference and surface tension
The cross-cell surface tension effect may be nullified or complemented by buoyancy, depending on whether the wetting or non-wetting fluid is denser. The critical value of 
$R$ corresponding to neutral stability is (see (4.1))
\begin{equation} R_c(k) = 1+\frac{k^2}{{\textit{Bo}}}+\frac{2 \alpha \cos \theta}{{\textit{Bo}} (b_0+ \alpha h)^2}. \end{equation}
The system is stable for all wavenumbers when 
$R < R_c(0)$, which is shown as a continuous blue line in figure 5. The results from § 3 for neutral stability for 
$k=0.2$ are plotted as black circles, showing good agreement. A comparison is also shown for 
$k=1.5$ as a broken blue line.
Neutral stability curves for (a) 
$\theta ={\rm \pi} /4$ and (b) 
$\theta =3{\rm \pi} /4$. Blue lines show the predictions of (4.8) for 
$k=0$ and 
$k=1.5$. Circles and crosses show the results from § 3 for 
$k=0.2$ and 
$k=1.5$, respectively. We use 
${\textit {Bo}}=5$, 
$h=0.5$, 
$b_0=0.3$, 
$M=2$, 
$\mathcal {A}=1$ and 
$C=1$.
Figure 6 shows the critical density ratio 
$R_c(0)$ as a function of the relative flux, and the channel inclination 
$\alpha$ in the case of a constant cell area, 
$0.2$. The interface height 
$h$ is obtained from (4.4). When 
$\alpha$ is small, the critical density ratio becomes independent of the relative flux (and hence the interface height 
$h$) because wetting effects become unimportant.
5. Conclusion
We have obtained the dispersion relation for the co-flow of two immiscible fluids in a Hele-Shaw cell with vertically varying gap width. The stability of the system is accurately predicted by the sign of the quantity
\begin{equation} J=R-1-{\textit{Bo}}^{{-}1}[k^2+2\alpha\cos \theta/b(h)^2]. \end{equation}
The last term, associated with channel wall inclination, represents the preference of the non-wetting fluid to occupy the thicker part of the channel. The interface is stable when the fluids have equal density and the wetting fluid occupies the thinner part of the channel. A density difference may complement or oppose this effect. We have obtained critical values of the density ratio, 
$R$, below which the system is stable to all wavenumbers. Our results also provide a basis for exploring the stability of important but more complex situations such as cells with elastic walls and cells whose vertical structure varies in the horizontal direction (e.g. Pihler-Puzović et al. Reference Pihler-Puzović, Périllat, Russell, Juel and Heil2013).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Coefficients for Frobenius’ method
In either fluid, the governing equation takes the form
\begin{align} &(\bar{b}^2 + a_4 \bar{b}^4 +a_6 \bar{b}^6) v_{\bar{b}\bar{b}} + (- \bar{b} + a_4 \bar{b}^3 +a_6 \bar{b}^5) v_{\bar{b}} +({-}3 + ({-}1-a_4) \bar{b}^2 \nonumber\\ &\quad -(a_4+a_6) \bar{b}^4 -a_6 \bar{b}^6) v = 0. \end{align}
The indicial polynomial is 
$n^2 - 2n -3 = 0$, which has solutions 
$n=3$ and 
$n=-1$. We write the first power series of 
$v(\bar {b})$ as
\begin{equation} \varPhi(x)=x^3 \sum_0^{\infty} P_{n} x^{n}, \end{equation}
with 
$P_0=1$ and the recurrence relation
\begin{equation} P_n \left[n^2+4n\right] + P_{n-2} \left[ n(n+2)a_4-1 \right] + P_{n-4}\left[a_6(n^2-2n) -a_4\right] -P_{n-6}a_6 = 0. \end{equation}The second independent power series solution is given by
\begin{equation} \varPsi(x) = \varPhi(x) \log x + x^{{-}1}\sum_0^{\infty} Q_{n} x^{n}, \end{equation}
where 
$Q_0= 16/(4a_4-1)$, 
$Q_2 = -Q_0/4$, 
$Q_4 =1$ and
\begin{align} 0 &= (n^2 -4n)Q_n +Q_{n-2}\left[a_4(n-2)(n-4)-1\right] + Q_{n-4}\left[a_6(n-4)(n-6)-a_4\right] \end{align}
\begin{align} & \quad- Q_{n-6} a_6+P_{n-4}(2n-4)+P_{n-6} a_4(2n-6) + P_{n-8}a_6(2n-10). \end{align}
Open access
