Bifurcations of drops and bubbles propagating in variable-depth Hele-Shaw channels

The steady propagation of air bubbles through a Hele-Shaw channel with either a rectangular or partially occluded cross section is known to exhibit solution multiplicity for steadily propagating bubbles, along with complicated transient behaviour where the bubble may visit several edge states or even change topology several times, before typically reaching its final propagation mode. Many of these phenomena can be observed both in experimental realisations and in numerical simulations based on simple Darcy models of flow and bubble propagation in a Hele-Shaw cell. In this paper, we investigate the corresponding problem for the propagation of a viscous drop (with viscosity ν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} relative to the surrounding fluid) using a Darcy model. We explore the effect of drop viscosity on the steady solution structure for drops in rectangular channels or with imposed height variations. Under the Darcy model in a uniform channel, steady solutions for bubbles map directly on to those for drops with any internal viscosity ν≠1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu \ne 1$$\end{document}. Hence, the solution multiplicity predicted for bubbles also occurs for drops, although for ν>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu >1$$\end{document}, the interface shape is reversed with inflection points appearing at the rear rather than the front of the drop. The equivalence between bubbles and drops breaks down for transient behaviour, at the introduction of any height variation, for multiple bodies of different viscosity ratios and for more detailed models which produce a more complicated flow in the interior of the drop. We show that the introduction of topography variations affects bubbles and drops differently, with very viscous drops preferentially moving towards more constricted regions of the channel. Both bubbles and drops can undergo transient behaviour which involves breakup into two almost equal bodies, which then symmetry break before either recombining or separating indefinitely.

bubbles, capsules and cells [3][4][5]. The small scales of these systems mean that flow is often at low Reynolds number, with nonlinearity arising primarily through the ability of two-fluid interfaces to deform. For strongly confined drops, the wetting properties of the fluids with the channel walls and any surfactants present may have a significant influence on the dynamics, along with surface tension acting on the deforming interface and viscous forces within each fluid domain. In many situations involving quasi two-dimensional drops, these ingredients combine to yield almost circular drops, which propagate steadily with little change in form. However, much more diverse behaviour can be observed in experiments where the channel topography is varied via a centred occlusion which penalises centred propagation (Fig. 1b), leading to a wide range of observable behaviour including symmetry breaking between steady states [6], oscillations relative to the moving tip of the finger or bubble [7][8][9], and bubble breakup followed by recoalescence or separation [10,11].
Full resolution of the dynamics of strongly confined bubbles or drops in 3D numerical simulations can be extremely challenging due to the disparate lengthscales in the film thickness, bubble diameters and channel dimensions. Most studies have focused on the equivalent of the experimentally observed solution branches: nearly circular drops, or smooth, half-width air fingers, each symmetric about the channel centre line and with no inflection points. Recent 3D studies for strongly confined drops [12,13] confirm that the films which confine the drops are not uniform in thickness. These thin films affect the drag on the drop [14][15][16][17] and have an important role in systems with surfactants as they dominate the surface area of the drop, and hence, any surface transport processes [18,19].
For air bubbles or fingers in channels with large aspect ratio, the Hele-Shaw approach is to replace the threedimensional free boundary value problem for fluid flow and interface motion with a reduced-dimensional description based on this large aspect ratio. In channels with constant rectangular cross sections and zero Reynolds number, the simplest reduced-dimensional model is the Darcy model comprising a two-dimensional, harmonically varying pressure in each fluid, or constant pressure inside any bubbles [1,20]. Inertia is neglected and the drop fills a fixed proportion of the channel height. The reduced model is second order in space, while the full Stokes equation is fourth order. As a result, not all boundary conditions can be applied. Typically (e.g. [21]) no-penetration conditions but not no-slip are applied at channel side walls, while at the fluid interface, a pressure jump is enforced relating to surface tension, which is a large aspect ratio equivalent of the normal stress condition. The tangential stress condition is neglected entirely and in the extension to the two-fluid case continuity of tangential velocity can also not be applied, thus, neglecting some important coupling for two-fluid or Marangoni cases [13,19,22]. In the absence of surface tension, [1,2] showed that steady propagation in the Darcy model yields a continuum of solutions, for both bubbles, with fixed finite volume, or for semi-infinite air fingers. The introduction of surface tension reduces this multiplicity via exponential asymptotics to a countable sequence of solution branches. For bubbles, the only stable solution is the single-tipped, symmetric solution branch calculated by [20]. The remaining branches become progressively more unstable and deformed, with each branch associated with an extra tip on the front of the bubble shape. Detailed 3D numerical calculations by [23] showed that this depth-averaged model accurately predicts the speed of the main solution branch for air fingers, provided that the channel aspect ratio W/H = α > 8.
Despite its simplicity, the Darcy model has proved remarkably accurate at predicting the diverse range of behaviour observed for air bubbles or fingers in channels with variable topography. This includes symmetry breaking for steady states, the onset and form of oscillations [21,24], and transient behaviour which involves bubble breakup and subsequent merging [8,10,11]. The model predicts each of these behaviours and also suggests that the complexity of the observed behaviour is particularly related to the first few, weakly unstable, solution branches in the Saffman-Taylor system. As the channel topograpy profile is varied, these solution branches can interact and create bistability between different steady propagation modes. Several unstable solution branches remain, and their unstable manifolds influence the transient behaviour under which bubbles can reach a stable, asymmetric solution branch via a process which surprisingly often involves topological changes. The subsequent two-bubble interaction is yet more sensitive to unstable edge states, with placement relative to these unstable states controlling whether bubbles will eventually recoalesce or separate indefinitely.
Based on the success of these models in predicting the dynamics of deformable bubbles, here, we explore how the richness of behaviour would be modified if the air phase is replaced with a second viscous fluid with dynamic viscosity ν relative to the surrounding, driving phase. The Darcy model fails to capture all the physical mechanisms of the system, with some of these deficiencies perhaps of greater relevance for two-fluid systems than for the case of air bubbles. However, this modelling approach has nonetheless proved remarkably effective at predicting the behaviour of bubbles and may still be of qualitative value in exploring viscous drop propagation. The key contributions of this paper concern how viscosity ratio qualitatively affects the steady and transient behaviour in channels with and without topography variations. We find that the full range of behaviour predicted by the Darcy model has a counterpart for viscous drops, and thus, the rich dynamical behaviour observed experimentally for bubbles is robust to drop viscosity within the validity of the model.
In Sect. 2, we summarise the depth-averaged model used for our multiphase system in the presence of variable topography and describe our finite-element implementation. In Sect. 3.1, we explore the steady solution structure for drops of viscosity ν relative to that of the surrounding fluid, in a channel of uniform height. In this case, we show that the solutions for drops and bubbles are equivalent, though appear at different flow rates for each viscosity. In Sect. 3.2, we explore the effect of centred expansions on the solution branches for steady propagation at each viscosity. For drops less viscous than the surrounding fluid (ν < 1), the behaviour is similar to that of air bubbles, with a central constriction leading to symmetry breaking of the main branch. The converse result appears for very viscous drops (ν > 1), with a centred constriction leading to stable symmetric propagation.
In Sect. 4, we explore transient behaviour from symmetric and asymmetric initial conditions. Again, we find the behaviour reverses when ν > 1, in that asymmetric initial conditions lead to asymmetric propagation modes for low-viscosity drops with central constrictions and high-viscosity drops with centred expansions. With appropriately selected initial conditions, we can also predict transient behaviour leading to breakup of the drop: a phenomenon often observed for experiments for bubbles. The details of the evolution towards breakup are different for air bubbles and viscous drops, with our model predicting that a viscous drop will create a thread connecting the two daughter drops which slowly thins while being drawn out by the bulk flow. Finally, in Sect. 5, we briefly simulate post-breakup behaviour and show that recoalescence and persistent separation are possible eventual outcomes.
We conclude with a critical discussion in Sect. 6 of how these behaviours link to the modelling assumptions, which phenomena could potentially be observed in multiphase experiments, and how depth-averaged models should be extended to better capture the full range of drop dynamics.

Leading order depth-averaged model in a channel of variable depth
We consider incompressible Stokes flow inside the region: where x * , y * and z * are dimensional coordinates as indicated in Fig. 1, W is the channel width and H its nominal height. The outer fluid has dynamic viscosity μ, while the drop has viscosity μν. The case ν = 0 corresponds to an air bubble, while we will refer to the cases 0 < ν < 1 and ν > 1 as low-and high-viscosity drops, respectively, relative to the viscosity of the outer phase.
We nondimensionalise x * and y * with respect to W/2 and z * with respect to H . We define α = W/H and operate under the assumption that α 1. The system is driven by an imposed pressure gradient −G * in the x-direction. We nondimensionalise p * with G * W/2 and the velocities with respect to the depth-averaged velocity In the Darcy model, the flow of each fluid is described by its leading order counterpart for large α. The pressure is then independent of z, and the depth-averaged velocity is u = −b 2 ∇ p in the outer fluid and u = −b 2 ∇ p/ν inside the drop. In each domain, mass conservation applies, so The flow is driven by the imposed pressure gradient: At the channel side walls, no-penetration conditions are applied, so p y = 0 at y = ±1.
In the limit of large α, pressure is the leading order contribution to the normal stress jump, and so the dynamic boundary condition can be replaced by the Young-Laplace equation, with a prescribed pressure jump which becomes where Q = μU * 0 /σ and σ is the coefficient of surface tension at the two-fluid interface. The right hand side of (1) represents the curvature of the interface, with contributions from the curvature of the drop as viewed from above, and also from the curvature of a semi-circular interface confined within the channel height. Note that this model does not include the small capillary number corrections derived by [25]. At leading order for Ca 1, this correction would have the effect of multiplying κ by π/4. This constant correction factor can be absorbed in our model by a rescaling of both Q and α and so does not affect the qualitative behaviour of either our steady or transient results.
If the drop fills the whole height of the channel, mass conservation implies that where U(t) = U (t)e x is the frame speed, e x is the unit vector in the x direction, n a unit normal vector pointing into the drop and R t is the velocity of a point on the drop interface in this moving frame. We can rearrange this expression so that n · ∇ p| fluid = − 1 b 2 (U · n + R t · n) at the interface and n · ∇ p| drop = − ν b 2 (U · n + R t · n) at the interface In our numerical method, we formulate the system in terms of the pressure, and hence, this boundary condition is the only equation where the viscosity ν arises. The channel topography function b is independent of x, and hence, the system can support steady propagation in the x direction. We take with s = 40 and w = 0.25. Solutions are shown in this paper for h = 0, h = 0.05 and h = −0.05, which correspond to a rectangular channel, or with a 5% centred constriction or expansion respectively. We take α = 40 throughout. No depth-averaged model will capture all the details of the full three-dimensional flow. This model neglects any films above and below the drop, and the low order (second-order two dimensional rather than biharmonic three dimensional) of the governing equations means that not all boundary conditions can be applied. We note in particular that continuity of tangential velocity cannot be applied at the interface, nor the balance of tangential stress.

Numerical methods
We use a finite-element method with a boundary-fitted mesh, similar to our previous work. The system is implemented in the finite-element library oomph-lib [26]. The drop and outer fluid domains are meshed using a triangular mesh with nodes coinciding at the drop boundary. Mesh deformation is handled by treating the mesh nodes as fixed points in a fictitious elastic solid, with stresses applied to this solid as a Lagrange multiplier to deform the interface as required.
For initial value problems and solutions for steady propagation, we work in a frame moving with the centre of mass of the drop. Hence, U(t) = U (t)e x , with U (t) varied dynamically so that the centroid position is constrained to always lie at x = 0. The domain is truncated at a distance L ahead of and behind the centroid, with p = 0 at x = L and p x = −1 at x = −L. We typically take L = 4, which is more than sufficient to contain the drop.
For two-fluid problems, there is no need to constrain the drop volume in time evolution. For steady calculations, we apply one extra constraint for each drop, 'trading' the mass conservation equation at a single node for the condition that the drop volume, calculated as follows: where S is the area of the drop as viewed in (x, y) space and ∂ S the drop interface, takes the required value.
Piecewise quadratic elements are used for the spatial discretisation, while unsteady calculations use the implicit BDF2 timestepper. Solutions are validated via convergence testing with respect to spatial mesh resolution, time stepping and domain truncation. Additional validation is provided by comparison to inviscid bubble solutions and via the mapping between solutions for drops of differing viscosity discussed in Sect. 3.1. We are able to solve , we see multiple 'tips', or nascent fingers, at the front of the bubble, with the branch number corresponding to the number of tips. These become more pronounced at lower Q. As ν increases towards ν = 1, these features become more pronounced at the same flow rate. For ν > 1, the situation is reversed. The tips appear at the rear of the drop and become less pronounced as ν increases for steady states directly, with initial guesses provided from the results of time-dependent simulations, from a variety of elliptical configurations or by parameter continuation. The system involves up to around 40,000 degrees of freedom, so a direct Newton solver is not prohibitive. For linear stability results, we formulate a generalised eigenvalue problem which is solved using the iterative Anasazi solver from Trilinos [27].

Solutions in rectangular channels
For the case of bubbles (ν = 0) steadily propagating in a rectangular channel, there is a countable sequence of solution branches. A representative sequence of drop shapes on these branches is shown in the top row of Fig. 2. All known solutions are symmetric about the channel centre line y = 0. The only stable solution branch corresponds to bubble shapes that are slightly pointed at the front and rounded at the rear, and elongated in the x-direction [28]. Subsequent branches become progressively more unstable and have increasingly distorted shapes. At low Q, these bubble shapes develop pronounced tips on the front of the bubble. Branch n (labelling the main stable branch as branch 1) develops n tips and is unstable to 2(n − 1) eigenmodes. All eigenmodes are either symmetric or anti-symmetric about y = 0, and each branch is unstable to n − 1 symmetric eigenmodes and n − 1 anti-symmetric eigenmodes. Similar solution multiplicity, involving a countable sequence of branches ( Fig. 2) with branch number corresponding to the number of tips (Fig. 3), exists for all non-unity viscosity ratios. The main qualitative difference is that the 'tips' appear on the rear of the drop when ν > 1. The same association of linear stability with tip number occurs for all viscosity ratios. When ν < 1, the stable branch corresponds to the fastest propagation, while for ν > 1, the stable branch is the slowest moving. In both cases, the stable shape is the narrowest, and U approaches 1 as the drop width or branch number increases.
Although the qualitative behaviour of the branches does not vary with viscosity, the choice of steady shapes at a fixed value of Q depend on the viscosity ratio. In fact under the Darcy model, we can find a non-trivial transformation between solutions for any non-unity viscosity ratio, but mapping to different Q as described below. As all solution branches appear to exist for all Q with no interactions between them, this mapping means that the same solution multiplicity occurs for each drop viscosity.
To derive this equivalency, we begin by considering the pressure inside the drop which obeys with its normal derivative on the drop boundary prescribed by the condition of steady propagation: For given U and ν, the internal pressure is uniquely specified by this boundary condition, up to an additive constant. Evidently, p = −νU x satisfies the required equations, and hence, we can write inside the droplet. The velocity inside the drop is simply u = U e x ; the constant p 0 does not affect the flow or drop shape. In the exterior of the drop,p = p + x satisfies From the kinematic condition, we have that The solution forp in the outer fluid is not trivial, but the problem is linear in U − 1, so we can let whereP is a function only of the interface shape andp 0 is a constant. All that remains is the dynamic boundary condition, which in the case b(y) = 1 reduces to Substituting the expressions for p outer and p inner from (2) and (3) into (4), we find Note that κ is also strictly a function of the interface shape. Variation around the interface in (5) arises only from x,P and κ. Hence, we can write (5) as follows: where F 0 may be a function of the global drop shape but is spatially constant, associated with addition of an arbitrary constant to the pressure jump p 0 −p 0 . The same drop shape can be a solution of the equations for different control parameters if U , ν and α 2 Q are varied in such a way to maintain the ratio between the three terms on the left-hand side of (6) which contain the only spatially varying terms in this equation.
As a result, an interface shape corresponding to a solution at ν 0 , U 0 , α 2 0 Q 0 is also a solution at a different viscosity ν, so long as U and α 2 Q are also varied to obey This result means that the same sequence of solution branches appears for each ν = 1 but are stretched in both Q and U (one example of this shown in Fig. 3). While respecting the restriction α 2 Q > 0, we may use this procedure to map solution branches with ν > 1 to each other and similarly for ν < 1. We can also map solutions with ν > 1 to solutions with ν < 1 and vice versa. This latter procedure generates a negative flow rate Q, but by reflecting the system in the line x = 0, we can recover a positive flow rate with the same U . This is why solutions at the same Q in Fig. 2 are reversed when ν > 1. We note that for each viscosity ratio, the only stable solutions in this figure are those in the far left column. Subsequent columns correspond to increasingly unstable solution branches which would not be expected to be observed in experiments.
In experiments, the prescribed flow rate Q is the control parameter, and U is a global solution measure. Under the transformation (7), solution branches at the same Q for one viscosity are mapped to different values of Q at a second viscosity. Hence, the selection of possible states in one experiment will depend on the viscosity.
For the case ν = 1 in the Darcy model, the background flow u = (1, 0) can be subtracted from both outer fluid and drop without changing the equations, if U → U − 1. Hence, steadily propagating states at ν = 1 with a background flow are equivalent to static states. These are circular drops with arbitrary offset, each with U = 1. No other simply connected steady or steadily propagating states exist in this case.
The equivalence between bubbles and drops in (7) relates to the solutions in the absence of surface tension given by Saffman and Taylor [2]. This result (7) is valid for multi-body steady states, so long as the two bodies have identical viscosity and surface tension. The equivalency (7) breaks down in several cases: unsteady motion, in the presence of height variation, and with more complex models such as the Brinkman equations where the reinstatement of continuity of tangential velocity and stress means that the flow inside the drop is not simply unidirectional. We will not consider these more complex models here. However, it seems likely that for sufficiently wide channels, the solution multiplicity predicted by the Darcy model will persist in some form.

Solutions in partially occluded channels
For bubbles, it is already known that the introduction of a centred occlusion influences the bifurcation structure [8,10,11,24,29]. Axially uniform height variations of as little as 2% are enough to induce pitchfork bifurcations and interactions between branches.
The topography profile b(y) features in the bulk equation and in the kinematic and dynamic boundary conditions. The dependence of the bulk and kinematic boundary conditions on b(y) relate to variations in viscous resistance in the channel and its effective permeability. The appearance of b(y) in the dynamic boundary condition represents a different physical mechanism reflecting the variations in the transverse curvature and hence surface tension of a semi-circular interface confined between the upper and lower plates.
The effect of b(y) on the viscous resistance breaks the equivalency between drops and bubbles. Specifically, with arbitrary b(y), the flow inside the drop is no longer unidirectional, and the excess pressure outside the drop is no longer proportional to U − 1. As a result, contributions appear in the dynamic boundary condition with up to five terms of different spatial dependence, too many to recover by rescaling of U and Q.
With the loss of this equivalency the introduction of a centred constriction affects bubbles and drops in different ways. For the case of bubbles (see e.g. [10]) or low-viscosity drops (the ν = 0.5 case in Fig. 4), the stable main branch interacts with the second solution branch. Symmetry breaking produces a stable asymmetric solution branch, while a limit point connects the next two unstable, symmetric solutions. As h is increased, the pitchfork bifurcation can move towards Q = 0, and for sufficiently high obstacles (e.g. for h = 0.05 in Fig. 4) the asymmetric, stable, solution branch has disconnected entirely from the symmetric branches at low Q which are now all unstable to symmetry-breaking perturbations.
For high-viscosity drops with ν > 1, the main, centred, solution branch exists for all Q with h = 0.05 with no loss of symmetry for this range of Q. In this case, the second and third solution branches interact at a limit point at Q ≈ 0.03. A supercritical pitchfork bifurcation occurs on the second solution branch before the limit point, producing a pair of asymmetric solution branches. However, as the second solution branch is itself doubly unstable, so also is the pair of asymmetric branches which emerge from the pitchfork.
This combination of a pitchfork bifurcation and a limit point, which occurs in both cases in Fig. 4 line y = 0. The main branch is linearly stable. The next branch is doubly unstable, with one symmetric and one anti-symmetric eigenmode. The third branch is quadruply unstable, with two symmetric and two anti-symmetric eigenmodes. It appears this sequence continues indefinitely, with each subsequent branch unstable to an extra symmetric and anti-symmetric eigenmode.
For two adjacent symmetric branches to meet at a limit point, one of the symmetric eigenmodes exchanges stability there. However, we must also exchange stability of one of the anti-symmetric eigenmodes. The simplest manifestation of this exchange of stability is via a pitchfork bifurcation, which must occur on one or other of the branches, before the limit point can occur. For the bubble-like ν = 0.5 case, the pitchfork emerges at very low h on  1 + ν), and the displacement in the x direction is recreated from observations of U the main branch, to produce a branch of asymmetric stable solutions which then persist for all flow rates. For the higher viscosity ν = 2 case, the pitchfork is associated with the interaction of two unstable branches and so does not produce a stable propagation mode.

A single drop or bubble
For arbitrary initial conditions, drops need not necessarily evolve directly to a stable state, especially when several steady solutions exist. In this section, we explore transient behaviour for drops at a range of viscosities for initial conditions corresponding to an offset, elongated drop (Fig. 5), and a centred, widened drop (Fig. 6).
For channels with a centred occlusion (h = 0.05), the offset low-viscosity drop (ν = 0.5) evolves to the stable asymmetric state, which is in fact very close to the initial condition. In contrast, for the higher viscosity case (ν = 2), the drop moves quickly to the centre of the channel (notably the rear of the drop moves over first) and adopts a stable propagation mode. Hence, it appears the low-viscosity drop seeks wider regions of the channel, whereas drops more viscous than the surrounding fluid prefer to occupy the more constricted regions.
We can confirm this preference by exploring channels with a centred expansion (h = −0.05), also shown in Fig.  5. In this case, the low-viscosity drop moves towards the centre of the channel, while high-viscosity drops migrate to an asymmetric steady state.
We also explore transient evolution from a symmetric initial state, where the possible migration to an asymmetric state is less straightforward. In Fig. 6, we choose an initial condition of a wide, centred ellipse, for a variety of viscosity ratios, and for constricted, rectangular and expanded channels. We show snapshots of the interface shape at regular intervals, with the displacement along the channel reconstructed from U (t). Drops either move towards a steady state, in which case the steady state itself is shown in the final column, or evolve towards breakup. Over For the rectangular channel (h = 0), the main branch solution is stable and symmetric for each drop viscosity. In all cases except ν = 4, the drop slowly narrows and appears to be evolving towards this stable solution. For ν = 4, the drop becomes horseshoe shaped, with growing fingers at its rear. The connecting region stretches and thins towards breakup.
With a centred constriction (h = 0.05), drops less viscous than the surrounding fluid develop growing tips on their front which eventually split the drop into two. Higher viscosity drops, with ν > 1, remain simply connected, elongating into the constricted central region. Finally, with a centred expansion (h = −0.05), the converse behaviour arises: low-viscosity drops elongate into the central region, while very viscous drops develop growing tips on their rear which will eventually split the drops into two.
For the case ν = 1, steadily propagating solutions are circular drops of arbitrary offset. Height variation breaks this degeneracy, and in fact, the shape of the drop is surprisingly strongly affected by the presence of height variation. The drop shape has two distinct tips and a pointed rear for h = 0.05, and a pointed front and flattened rear for h = −0.05.

Drop breakup
The breakup of a drop into two equal parts is observed for bubbles [8,10,11], and for both low-and high-viscosity drops (Fig. 6). Evolution towards breakup is one case where we might expect different behaviour for bubbles with ν = 0 to very low-viscosity drops. If ν = 0, the bubble simply acts as a volume constraint, with no internal flow. The depth-averaged interface, hence, self-intersects in finite time, at which point breakup is inevitable.
Simulations for drops of any viscosity also show evolution towards breakup via the same symmetric tip-splitting as observed for bubbles. However, when the drop has internal viscosity, the breakup does not occur immediately. Instead, a persistent strand connects the two daughter drops. This strand is stretched by the flow and stretched considerably. The strand thins during this process and in some cases appears to develop beads (e.g. for the ν = 4, h = 0 case in Fig. 6). The meshing demands and computional complexity increase rapidly as the strand thins, and it is not clear whether a simply connected solution exists to the model equations for all times. The last state computed for the simply connected drops is shown in Fig. 7.

Multiple drops
The finite-element discretisation we are using for our system does not naturally allow for topographical changes. Hence, changes in the drop topography require specific intervention: we restart the simulation with a new initial condition corresponding to the post-breakup state. As the flow is at zero Reynolds number, the only initial condition needed is a specification of the interface shape, obtained from the last simply connected state by removing the thread at the point where it joins each drop. The last simply connected single-drop state is shown in blue, and the two-drop state shown in bold, in Fig. 7. In these cases, the one-drop simulation began from a centred initial condition, and the subsequent breakup yields an almost symmetric two-drop state, with one daughter drop on each side of the channel. In both experiments and numerical simulations, this post-breakup condition will never be perfectly symmetric, whether due to mesh-induced asymmetries or unavoidable experimental imperfections.
Just after breakup, the drop shape is sharply pointed in the region near where the thread was removed. This sharp region quickly smooths with little effect on the rest of the drop shape, and the pair of drops appears to be attracted to a symmetric two-drop steady state with one drop on each side of the channel. However, this two-drop steady state is unstable, and in each case, we observe that one drop moves ahead of the other, breaking the symmetry of the steady state. The drops are somewhat attracted to each other in both cases, but for these particular examples, different outcomes arise from the two cases. The low-viscosity drops eventually separate, while the high viscosity drops merge.
For air bubbles, two-bubble steady states have a rich bifurcation structure of their own [11], with unstable twobubble states act as edge states in organising the transient dynamics and ultimately determining whether bubbles merge or separate. Depending on the system parameters and the size ratio post-breakup, separation or recoalescence can be observed in experiments and predicted in Darcy-based models, with good qualitative agreement between the two.
For drops of equal size and viscosity in a rectangular channel, a two-drop symmetric steady state can be constructed by placing two drops in two half-width channels (first column in Fig. 8). This construction follows from the Darcy model as the channel side walls at y = ±1 and channel centre line y = 0 are each simply lines where p y = 0. These two-drop states are unstable to symmetry-breaking perturbations, with the eigenmode corresponding to one drop moving ahead and the other behind.
A branch of two-drop states can be tracked by continuation from the two-equal drop case with h = 0. Continuation in h to ±0.05 results in only very minor changes to the drop shape (first and second columns in Fig. 8). If the drops are identical in volume, viscosity and surface tension, then the resulting steady state remains symmetric, but is unstable to symmetry-breaking perturbations.
We can use continuation methods from these known two-drop states to map out solution branches for the steady propagation of pairs of non-identical drops, for example, by varying the viscosity of one drop while holding the other fixed (subsequent columns in Fig. 8). The two-drop state at h = 0 is unstable, and unless any bifurcations occur, this stability property will still hold for states obtained from simple continuation. The propagation speed of a single drop reduces with viscosity (in this model, a circular drop in an unbounded cell moves with U = 2/(1 + ν)). The propagation speed of a single drop also depends on size and flow rate, though more weakly than on viscosity. This speed difference means that steadily propagation collections of drops with different viscosities rely on strong interactions between the drops. We can construct steady states with drops of slightly different viscosity, or varying surface tension or size. The two-equal-drop steadily propagating state can only bear a minor viscosity difference (see sequence in each row of Fig. 8), as the interaction and deformation of the drops is required to counteract a substantial difference in the pressure gradients required to drive the same internal flow, and hence, the same propagation speed as viscosity varies.

Discussion
In this paper, we investigated the dynamics of drops of arbitrary viscosity in a Hele-Shaw channel with and without imposed axially uniform variations in the channel topography. Our analysis was based on a simple depth-averaged Darcy model which has previously been shown to qualitatively capture many phenomena observed in experiments involving the propagation of inviscid air bubbles in the same system. For the Darcy model, steadily propagating bubbles and droplets of any viscosity ratio ν = 1 are equivalent under a transformation involving the independent flow rate Q and the solution measure U . Physically, this is because the model predicts unidirectional flow within any steady propagating drop, which matches the driving pressure gradient. This symmetry is also responsible for features such as tips appearing on the front of low viscosity drops but the rear of higher viscosity drops. As a result, when h = 0, the same solution structure of disconnected branches exists for each viscosity ratio, with the exception of ν = 1. This equivalence between drops and bubbles is broken by transient behaviour, variations in the channel topography, multiple drops with different viscosity or surface tension, or an alternative flow description which predict a non-unidirectional flow within a steadily propagating drop.
We focused mainly on the effect of varying channel topography, in the form of a centred, axially uniform, constriction or expansion in the central part of the channel. These height variations affect bubbles and drops in different ways: with central constrictions, low viscosity drops (ν < 1) more readily symmetry break, and in fact, a stable asymmetric solution branches exist for high obstacles at all flow rates. High-viscosity drops (ν > 1) migrate towards constricted parts of the channel, where the flow is generally slower. In this case, the stable branch is symmetric and larger drops elongate along this constricted region if necessary. The converse results are obtained for centred expansions: less viscous drops prefer to occupy the centred expansion region, while very viscous drops move away from the centre of the channel.
For low-viscosity drops in channels with centred constrictions, or high-viscosity drops in the case of centred expansions, an ubiquitous phenomenon observed in numerical experiments is the breakup of a large drop into smaller ones, driven by tip splitting at the front or rear for low-and high-viscosity drops, respectively. The resulting two-drop state is symmetric about the channel centre line but unstable to symmetry breaking perturbations.
For air bubbles, the post-breakup dynamics is sensitive to unstable two-bubble steady states which act as edge states in organising trajectories within a high-dimensional phase space and hence determining whether bubbles eventually merge or separate indefinitely, with similar dynamics occurring in both experiments and in numerical simulations based on the Darcy model [8,11]. We confirmed that both of these ultimate outcomes are predicted to arise for viscous drops. If drops interact sufficiently strongly that one crosses the central part of the channel, both drops reach the same side of the channel and may then move towards each other until coalescence would occur. Alternatively, with weaker interactions, the drops remain on opposite sides of the channel and slowly drift apart.
The Darcy model is of lower order than the underlying Stokes equations, and hence, the boundary conditions corresponding to continuity of tangential velocity and stress cannot be applied. It is possible to incorporate these conditions via an improved depth-averaged framework similar to the Brinkman equations [22]. Our preliminary results for propagation of drops in an unbounded Hele-Shaw cell with constant height suggest that at least some of the multiple solution branches still exist in this model. An exploration of the bifurcation structure and the effect of topography on drops and bubbles as predicted by this modelling approach is being planned.
The approach to breakup or coalescence is a demanding test for any depth-averaged model: the separation of scales which describes a single undeformed drop-that flow variations are rapid in the z-direction, and of the same slow order in the x-and y-directions-breaks down, so should be approached in a different limit [30]. In our model, we assumed that the drop fills the whole height of the channel. This assumption prevents the Rayleigh-Plateau instability, which would be expected to lead to faster pinchoff of thin viscous threads and also means that the exterior fluid cannot pass above and below the thread and must instead drag the thread with the bulk flow. The evolution of this viscous thread is one of the main differences between low-viscosity drops and strictly inviscid bubbles. Bubbles break in finite time as the mesh self-intersects. The approach to coalescence for bubbles is slower, and in our simulations [11], coalescence is modelled by a simple distance tolerance. This distance criterion will not capture the details of coalescence itself, but mimics to some extent the Rayleigh-Plateau instability and does appear to provide reasonable accuracy for coupling to the macroscale description of the system.
Experimental implementation of the two-fluid system described here is not necessarily a simple extension of the air bubble setup used in studies such as [8,10,11,24]. In particular, we have assumed that the surrounding fluid perfectly wets the channel walls and hence that there is no contact line, and our model does not account for surfactants. It may be possible to include the effects of surfactant transport, or of generally less-passive confining films, by developing lubrication equations to describe the evolution of a film height h f (x, y, t), which could be coupled to a depth-averaged description for the flow inside each drop. However, the required asymptotic matching at the edge of the drop to the bulk flow in the outer fluid would be non-trivial, and the predictions of such a model may well be sensitive to the assumptions made in this region.
One final consequence of equivalence between bubbles and drops in the Darcy model is that it cannot produce steady states for the propagation of double-emulsion drops. As described in Sect. 3.1, the flow inside the outer drop in a steadily propagating state is undirectional and the pressure is linear in x. An inner drop feels only this linear pressure gradient and no effect of the shape of the outer drop. Steady states can be created only if the viscosities of inner and outer drops match, but even then the system is degenerated as the inner drop could be placed anywhere inside the outer drop. Double emulsion drops are often used in microfluidic applications and evidently can propagate steadily. These double-emulsion drops are of course three dimensional and usually involve surfactants. However, even a relatively crude depth-averaged model may prove useful, for example, in exploring the physical mechanisms underlying the non-trivial interactions observed experimentally in trains of double-emulsion drops [31]. exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.