# Spreading or contraction of viscous drops between plates: single, multiple
or annular drops

H. K. Moffatt1 Howard Guest2 Herbert E. Huppert3 1Department of Applied
Mathematics and Theoretical Physics,
Wilberforce Road, Cambridge CB3 0WA, UK
2 St John’s Innovation Centre, Cowley Rd, Milton, Cambridge CB4 0WS
3 Institute of Theoretical Geophysics, King’s College, Cambridge CB2 1ST

###### Abstract

The behaviour of a viscous drop squeezed between two horizontal planes (a
contracting Hele-Shaw cell) is treated by both theory and and experiment. When
the squeezing force $F$ is constant and surface tension is neglected, the
theory predicts ultimate growth of the radius $a\sim t^{1/8}$ with time $t$.
This theory is first reviewed and found to be in excellent agreement with
experiment. Surface tension at the drop boundary reduces the interior
pressure, and this effect is included in the analysis, although it is
negligibly small in the squeezing experiments. An initially elliptic drop
tends to become circular as $t$ increases. More generally, the circular
evolution is found to be stable under small perturbations. If, on the other
hand, the force is reversed ($F<0$), so that the plates are drawn apart (the
‘contraction problem’), the boundary of the drop is subject to a fingering
instability on a scale determined by surface tension. The effect of a trapped
air bubble at the centre of the drop is then considered. The annular evolution
of the drop under constant squeezing is still found to follow a ‘one-eighth’
power law, but this is unstable, the instability originating at the boundary
of the air bubble, i.e. the inner boundary of the annulus. The air bubble is
realised experimentally in two ways: first by simply starting with the drop in
the form of an annulus, as nearly circular as possible; and second by forcing
four initially separate drops to expand and merge, a process that involves the
resolution of ‘contact singularities’ by surface tension. If the plates are
drawn apart, the evolution is still subject to the fingering instability
driven from the outer boundary of the annulus. This instability is realised
experimentally by levering the plates apart at one corner: fingering develops
at the outer boundary and spreads rapidly to the interior as the levering is
slowly increased. At a later stage, before ultimate rupture of the film and
complete separation of the plates, fingering spreads also from the boundary of
any interior trapped air bubble, and small cavitation bubbles appear in the
very low pressure region far from the point of leverage. This exotic behaviour
is discussed in the light of the foregoing theoretical analysis.

Keywords: Hele-Shaw cell, one-eighth power law, elliptic drop, adhesion,
lubrication, finite-time singularity, fingering instability, annular
evolution, contact singularity, cavitation

##  1 Introduction

The behaviour of a drop or film of viscous liquid trapped between two
horizontal plates (a Hele-Shaw cell), when the plates are squeezed together by
a force $F(t)$, presents a problem of great practical importance in the
adhesion industry. In this context, investigations date back to Stefan (1874),
who derived the relation between $F(t)$ and the separation of the plates
$h(t)$ when the upper plate is a circular disc that is impelled downwards
towards a fixed substrate, and the film extends to the disc boundary; the same
situation was treated in the more accessible review by Bikerman (1947).

When the drop does not extend to the plate boundaries, and when the squeezing
force $F$ is positive and constant (as is the case in figure 1), it is known
on the basis of lubrication theory that, when surface tension is neglected,
the separation $h(t)\sim t^{-1/4}$ (Moffatt 1977), and the corresponding drop
radius $a(t)$ increases like $a(t)\sim t^{1/8}$ (while the gap contracts, the
drop radius expands). This ‘one-eighth’ power law also applies to the
gravitational spread of a viscous liquid on a horizontal plane due to the
gradient of the free surface (Huppert 1982), with important implications for
the spread of lava from volcanic eruptions (Huppert et al. 1982). For the
squeezing problem, the effect of surface tension, which provides a pressure
jump across the air-liquid interface, has been considered by Ward (2006) (see
§2.1 below, and also §4 where stability of the basic state is considered).
Surface tension is of critical importance when the drop is considered as
forming a ‘liquid bridge’ between the two plates (see, for example, Gaudet,
McKinley & Stone 1996).

(a) $\tau=60$ (b) $\tau=172$ (c) $\tau=321$

Figure 1: Photographs of an expanding viscous drop placed between two
horizontal glass plates, the upper plate being allowed to descend under its
own weight; the drop is illuminated from below and viewed from above, and the
shots were taken at the dimensionless times $\tau=t/t_{0}$ indicated [here
$t_{0}=2.15$s, defined by (11)]. The background grid (which may be seen by
zooming to expand the figure) covers a 130 mm square in each panel, and allows
measurement of the radius $a(\tau)$ to within $\pm 0.1$mm. Note how the
opacity decreases as the layer thickness $h(\tau)$ decreases. Figure 2:
Continuation of figure 1, with levering at lower right-hand corner, in order
to gradually open the gap between the plates; the leverage here was increased
slowly for about one minute, before ultimate rupture of the film; (a) early
stage of fingering instability which starts from the ‘south-east’ sector of
the boundary; (b) tip-splitting and side-branching begins as fingers grow;
(c-e) the process continues, and in (e) cavitation bubbles can be detected in
the ‘north-west’ area of the drop; (f) fingering fills the fluid domain, and a
‘ridge’ is evident where fingers emanating from the north-west sector of the
boundary impact those from the south-east.

When the force $F$ is negative, i.e. when the plates are pulled apart with
initial conditions $a(0)=a_{0},\,h(0)=h_{0}\ll a_{0}$, then (see §2.2)
lubrication theory leads to the results $h(t)=h_{0}(1-t/|t_{0}|)^{-1/4}$ and
$a(t)=a_{0}(1-t/|t_{0}|)^{1/8}$, where $|t_{0}|$ is a characteristic time
proportional to $|F|/\mu$ and $\mu$ is the dynamic viscosity of the fluid.
This indicates an incipient ‘finite-time singularity’ at time $t=|t_{0}|$
when, in effect, the plates separate completely; but lubrication theory ceases
to be valid when the critical time is approached, specifically when
$1-t/|t_{0}|\sim(h_{0}/a_{0})^{8/3}$.

There are additional complications when $F\\!<\\!0$, which have been reviewed
by Gay (2002). First, there is an instability due to the suction of air into
the region occupied by the more viscous liquid as its area decreases (see
figure 2, which shows the effect of levering the plates apart at one corner).
This is a fingering instability, as first described by Saffman & Taylor
(1958), and as encountered in the Hele-Shaw experiments of Rabaut et al.
(1988), different only in that in these experiments air was injected at high
pressure into the lower-pressure viscous layer. Following the computational
investigation of Kelly & Hinch (1997), it was shown by Tanveer (2000) that the
stability problem is ill-posed in this sort of situation, but that the
potential development of singularities can be controlled by surface tension;
this is consistent with the analysis of the fingering instability that we
present in §5. Fractal patterns generated by fingering in a separating-plate
experiment have been compared with the results of statistical simulations by
La Roche et al. (1991). We note that fingering can be controlled to some
extent, not only by surface tension, but also through replacement of the
plates by elastic membranes, a situation of potential importance in biological
contexts (Pihler-Puzović et al. 2012).

Secondly, cavitation bubbles may appear in the course of plate separation,
presumably wherever the pressure would otherwise fall below the vapour
pressure of the viscous liquid; cavitation bubbles may be detected in figure
2(e). a phenomenon studied in detail by Lakrout et al. (1999) using a ‘probe-
tack apparatus’ and acrylic pressure-sensitive adhesives. When varying the
traction velocity, there is a transition from a fingering regime to a
cavitation regime, which has been investigated by Poivet et al. (2004). Again,
the effect of surface tension, which can counteract the traction force when
the gap is sufficiently small and is therefore crucial for effective adhesion,
has been treated by Ward (2011). This also will be reviewed in §2.2, and again
in the stability analysis of §5.

Figure 3: Evolution of an annular drop; stills at time $t$ (seconds) as
indicated for panels (a-f). Panel (a): the initial, nearly circular, annular
drop; (b) eruption of a secondary bubble from where the internal curvature is
maximal; the outer boundary is only weakly perturbed; (c) tip-splitting of the
secondary bubble, which has just penetrated the outer boundary; note that the
colour shows that a very thin ‘wetting’ layer of treacle has been left on both
plates where the secondary bubble has evolved; (d) eruption of the secondary
bubble to the exterior as it detaches from the primary bubble; (e) eruption
nearly complete; (f) eruption of tertiary bubbles from the primary bubble; (g)
levering at lower right-hand corner of the plates is introduced, showing early
development of fingering instability from the outer boundary; (h) the fingers
advance showing significant tip-splitting and side-branching; puncturing is
imminent; (i) a finger punctures the tertiary bubble (at the point marked by
the arrow) allowing equalisation of pressure with the exterior; (j) no further
puncturing occurs, but a ridge is established, separating the ‘invading’ and
‘defending’ cohorts of fingers; (k) fingering spreads right round the outer
boundary and cavitation bubbles appear ahead of the advancing front; (l) the
fingering pattern is now fully established, and the ridge remains prominent,
with a small gap only at the puncture point. Figure 4: Expansion and merging
of four drops of treacle between two horizontal plates; time $t$ in seconds as
indicated. An air bubble is trapped, but does not relax to circular form; the
volume of air in the bubble (or bubbles) is constant for so long as they
remain trapped. Panel (a): initial configuration; (b) the four drops have
expanded, and a first contact is imminent; (c) first two contacts have
occurred almost simultaneously, forming cusps which immediately resolve
themselves through the effect of surface tension; a third contact is imminent;
(d) the air bubble is now trapped, and its volume remains constant; (e) the
fourth cusp is resolved by air pushing into the surrounding treacle forming a
secondary bubble, while the combined volume of primary and secondary bubbles
remains constant; (f) the secondary bubble expands and (g) penetrates the
outer boundary; (h) the secondary bubble drains to the exterior, and a
tertiary bubble begins to form from the residual cusp on the primary bubble;
(i) the secondary bubble has escaped and the tertiary bubble grows very
slowly.  Figure 5: Continuation of figure 4. Panels (a-c): the tertiary bubble
grows extremely slowly for over three hours and appears to attain a quasi-
equilibrium. Panels (d-i): development of fingering instability induced by
levering the plates apart with a knife inserted at the lower right-hand
corner; the levering was gradually increased for about one minute before final
rupture; (d) early stage of fingering instability at the nearly circular outer
boundary of the drop; (e) early stage of tip-splitting; (f) early stage of
side-branching and growth of ‘defending’ fingers from the bubble boundary; (g)
advanced stage of fingering while the defending fingers retreat and one of the
‘invading’ fingers punctures a defending finger providing a pressure-
equalising path from the bubble to the exterior; (g) invading fingers advance
round the boundary, and cavitation bubbles appear ahead of the advancing
front; (i) ultimate stage of fingering, when the treacle is confined to the
narrow tree-branches and twigs that now separate the fingers. Note here the
near-pentagonal ridge that separates the invading fingers from the defending
fingers, with again a single breach in the south-east sector of the ridge.

There have been a number of investigations of the behaviour when the initial
shape of the drop is non-circular. This problem has been previously treated by
Savina & Nepomnyashchy (2015) (see also Savina 2021) by means of conformal
mapping, and using the Schwarz function (Davis 1974, Gustafsson & Vasiliev
2006) as earlier proposed for the two-phase (Muskat) problem by Howison
(2000). In this way, they found a solution for the case of an initially
elliptic drop (see §4.1 for a simpler treatment of this particular problem).
Savina & Nepomnyashchy (2015) also treated the problem of a nearly circular
hypotrochoid of three-fold symmetry, and came to the surprising conclusion
that, under squeezing and neglecting surface tension, the drop shape ‘evolves
towards a deltoid’, a curve with three cusps. This goes against physical
intuition, and our stability analysis in §4.2 comes to the opposite
conclusion, namely, that quite generally, under squeezing with a constant
force, any perturbation from the circular shape is eliminated as the area of
the drop increases. This behaviour is consistent with the experiments
described in §7.

When an air bubble is trapped, either accidentally or deliberately, inside the
viscous drop, new problems arise. When the viscous liquid has the form of an
annulus trapping a circular air bubble, the ‘basic-state’ evolution can be
easily determined and still follows a one-eighth power law evolution (§3); but
we shall find that this basic state, even under squeezing, is unstable, the
instability being driven (in a sense that will be quantified) from the bubble
boundary (§6). For the case of contraction ($F\\!<\\!0$), we shall find that
the fingering instability persists and is driven primarily from the outer
boundary.

The details of the experiments providing the photographs of Figures 1-5 will
be described in §7. These home experiments, conducted under lockdown
conditions, replicate lecture-room demonstrations that we have used over many
years to illustrate the one-eighth expansion law when $F\\!>\\!0$ and the
fingering instability when $F\\!<\\!0$. We used Lyle’s black treacle as the
working fluid, which allows for good colour contrast. When a single drop of
this fluid is placed on a fixed horizontal glass plate, and a second glass
plate is placed on the drop and allowed to descend under its own weight, the
drop expands and soon adopts an accurately circular form; the slow expansion
can be allowed to continue for several hours. The situation $F\\!<\\!0$ can
then be achieved by gently levering the plates apart at one corner. The
separation $h({\bf x},t)$ then varies linearly in ${\bf x}=(x,y)$, but the
fingering instability, non-uniform around the drop boundary, is nevertheless
well illustrated (figure 2). Furthermore, cavitation bubbles can be observed
where the pressure is minimal, and explained at least qualitatively in terms
of the ‘hinged plate’ problem described in §2 of Moffatt (1964).

The trapping of a large air bubble was achieved experimentally in two quite
different ways, first by starting with a drop in the form of a circular
annulus (see figure 3), and second, by starting with four disjoint drops
(figures 4 and 5), and forcing these to expand under the weight of the upper
plate so that they merge, trapping an air bubble that is very non-circular
from the outset. In both cases, the instability that is predicted
theoretically manifests itself through the ejection of a secondary bubble from
the primary bubble, allowing some air to escape. The subsequent evolution is
extremely slow, but the boundary of the air bubble remains very non-circular
throughout. When the plates are levered apart from one corner, the fingering
instability appears first from the external boundary, and later from the
internal bubble boundary, and a ridge develops, separating these two fingering
‘cohorts’.

##  2 The basic state

###  2.1 The squeezing problem $(F\\!>\\!0)$

Consider what happens when a drop of viscous liquid of volume $V$ and dynamic
viscosity $\mu$ is placed between two horizontal plates $z=0$ and $z=h(t)>0$,
the upper plate being impelled towards the fixed lower plate by a constant
force $F$. It is observed in this situation that, when viewed from above, the
drop tends to become of circular form of expanding radius $r=a(t)$, where,
with $h(0)=h_{0},a(0)=a_{0}$, by conservation of volume,

$\pi\,a^{2}(t)\,h(t)=\pi\,a_{0}^{2}\,h_{0}=V=\textnormal{const.}$ (1)

We shall describe this evolving circular situation, shown in figure 1, as the
basic state. In this state, provided $h(t)\\!\ll\\!a(t)$, standard lubrication
theory implies that (i) the pressure $p(r,t)$ is independent of $z$ and (ii)
the radial velocity component $u(r,z,t)$ satisfies $|\partial u/\partial
z|\gg|\partial u/\partial r|$ and is driven by the pressure gradient according
to the equation

$\frac{\partial p}{\partial r}=\mu\frac{\partial^{2}u}{\partial\,z\,^{2}}\,,$
(2)

with boundary conditions $u(r,0,t)=u$($r,h(t),t$)$=0$. Hence

$u(r,z,t)=\frac{1}{2\mu}\frac{\partial p}{\partial r}\,z\,(z-h)\,,$ (3)

and, averaging over $z$,

$\bar{u}(r,t)\equiv\frac{1}{h}\int_{0}^{h}u(r,z,t)\,\textnormal{d}z=-\frac{h^{2}}{12\mu}\,\frac{\partial
p}{\partial r}\,.$ (4)

Mass conservation and incompressibility indicate that the radial flux
$Q(r,t)=2\pi r\,h\,\bar{u}(r,t)$ at each radius $r$ is equal to the flux $-\pi
r^{2}\,\textnormal{d}h/\textnormal{d}t$ driven downwards by the upper plate.
It follows that

$\frac{\partial p}{\partial
r}=\frac{6\,\mu\,r}{h^{3}}\,\frac{\textnormal{d}h}{\textnormal{d}t}$ (5)

and we note that

$\nabla^{2}p=\frac{1}{r}\frac{\partial}{\partial r}r\,\frac{\partial
p}{\partial
r}=\frac{12\mu}{h^{3}}\,\frac{\textnormal{d}h}{\textnormal{d}t}\,.$ (6)

This is the particular form of ‘Reynolds’ equation’ applicable to this problem
(see for example §3.6 of Moffatt 1977).

Let $\gamma$ denote the surface tension at the liquid/air boundary. Assuming a
circular meniscus with zero contact angle, there is then a jump of pressure
$\gamma/h(t)$ with corresponding suction across this boundary. The boundary
condition is then

$p=p_{a}-\gamma/h(t)\quad\textnormal{at}\quad r=a(t)\,,$ (7)

where $p_{a}$ is the atmospheric pressure. With this boundary condition the
pressure in the drop is then given by

$p(r,t)-p_{a}=\frac{3\mu}{h^{3}(t)}\frac{\textnormal{d}h}{\textnormal{d}t}\left(r^{2}-a^{2}(t)\right)-\frac{\gamma}{h(t)}=p_{0}(r,t),\,\,\textnormal{say}\,.$
(8)

The total vertical force exerted by the fluid on the upper plate is then

$F\\!_{p}=2\pi\int_{0}^{a}p_{0}(r,t)\,r\,\textnormal{d}r=-\frac{3\pi\,\mu}{2h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\,a^{4}-\frac{\pi\gamma}{h(t)}a^{2}\,.$
(9)

Using (1) and neglecting plate inertia, the force balance $F=F\\!_{p}$, we
then find

$F=\frac{3\mu
V^{2}}{8\pi}\frac{\textnormal{d}}{\textnormal{d}t}\left(\frac{1}{h^{4}}\right)-\frac{\gamma
V}{h^{2}}\,.$ (10)

With dimensionless variables

$X=h_{0}^{2}/h^{2}(t)=a^{4}(t)/a_{0}^{4}\quad\,\,\textnormal{and}\,\,\
\tau=t/t_{0}\,,\quad\textnormal{where}\quad
t_{0}=\frac{3}{8\pi}\frac{\mu}{F}\,\frac{V^{2}}{h_{0}^{4}}=\frac{3\pi^{3}}{8}\frac{\mu}{F}\frac{a_{0}^{8}}{V^{2}}\,,$
(11)

equation (10) takes the convenient dimensionless form

$\textnormal{d}X^{2}/\textnormal{d}\tau\equiv
2X\,\textnormal{d}X/\textnormal{d}\tau=1+\lambda\,X\quad\textnormal{with}\quad
X(0)=1,$ (12)

where

$\lambda=(\gamma/F)\left(V/h_{0}^{2}\right)=(\gamma/F)\left(\pi^{2}a_{0}^{4}/V\right)\,.$
(13)

Here $\lambda$ is the relevant dimensionless measure of surface tension. It is
small $(\sim 0.003)$ in the experiments described in §7, but surface tension
effects can evidently still become important as and when $X$ increases to
O$\left(\lambda^{-1}\right)$, i.e. when $a/a_{0}$ increases to
O$\left(\lambda^{-1/4}\right)$.

(a) (b)

Figure 6: (a) Plot of $X^{2}(\tau;\lambda)=a^{8}(t)/a_{0}^{8}$, with
$\tau\equiv t/t_{0}$, as determined by (11) and (14), for
$\lambda=0.001,0.01,0.02$; note that $X^{2}\sim 1+\tau$ as $\lambda\rightarrow
0$; (b) plot of $X^{2}(\tau;|\lambda|)$, with $\tau=t/|t_{0}|$, for the case
of contraction $F<0$ [as given by (19)]; here $X^{2}\sim 1-\tau$ as
$\lambda\rightarrow 0$,

(a) (b)

Figure 7:  (a) Plot of $\hat{a}(\tau)\equiv a(t)/a_{0}$ for a squeezed drop;
the red curve shows the function defined by (16); the dotted lines correspond
to the three images shown in figure 1; (b) plot of $\hat{a}^{8}(\tau)$,
showing clearly the one-eighth power-law over the full range of $\tau$.

Equation (12) integrates to give

$\tau=\frac{2}{\lambda^{2}}\left((X-1)\lambda-\log{\frac{1+\lambda
X}{1+\lambda}}\right)\,,$ (14)

thus determining $X$ implicitly as a function of $\tau$ and $\lambda$ (cf.
Ward 2006). For $\lambda\ll 1$, this gives

$\tau=X^{2}-1-\textstyle{\frac{2}{3}}\displaystyle\left(X^{3}-1\right)\lambda+\textnormal{O}\left(\lambda^{2}\right)\,.$
(15)

Figure 6(a) shows $X^{2}(\tau)$ as determined by (14) for three small values
of $\lambda$; as expected from (15), the curves all asymptote to $X^{2}\sim
1+\tau$ as $\tau\rightarrow 0$.

If surface tension is neglected (i.e. $\lambda=0$), then (15) gives
$X(\tau)=(1+\tau)^{1/2}$ or equivalently

$h(\tau)=h_{0}(1+\tau)^{-1/4}\,,\quad a(\tau)=a_{0}(1+\tau)^{1/8}\,.$ (16)

The result for $a(\tau)$ is the ‘one-eighth power law’ (as indicated in the
introduction). Figures 1 and 7 show the results of a simple experiment
(details in §7); the one-eighth power-law shows up quite accurately in this
experiment, for which $t_{0}\approx 2.15$s and $\lambda\approx 3\times
10^{-3}$. The time taken to approach this one-eighth power-law is evidently
very short in this experiment. This ’adjustment time’ has been recently
investigated by Ball & Huppert (2019) and Webber & Huppert (2020) for the case
of gravitational spread of a viscous fluid (Huppert 1982). They have found
that this time, somewhat independent of initial conditions, can be just
fractions of a second in laboratory experiments, but many days in geological
situations, where the extremely viscous lava domes from volcanic eruptions can
eventually spread horizontally to hundreds of kilometres.

###  2.2 The contraction problem $F<0$

If $F<0$ (i.e. if the upper plate is pulled away from the lower plate with a
constant traction force $|F|$) leading to contraction of the area of the drop,
then $t_{0}<0$, and it is natural therefore to redefine $\tau$ as
$\tau=t/|t_{0}|$; with this modification, (12) becomes

$\textnormal{d}X^{2}/\textnormal{d}\tau\equiv
2X\,\textnormal{d}X/\textnormal{d}\tau=-1+|\lambda|\,X\quad\textnormal{with
again}\quad X(0)=1.$ (17)

Now a steady state with $X=1$ is possible if $|\lambda|=1$; but, if
$|\lambda|\neq 1$,

$\left|\textnormal{d}X^{2}/\textnormal{d}\tau\right|_{\tau=0}=-1+|\lambda|\,,$
(18)

so that $X^{2}$ increases or decreases from 1 according as $|\lambda|\gtrless
1$. If $|\lambda|>1$, the reduction in pressure in the drop due to surface
tension is sufficient to overcome the traction force, thus decreasing
$h(\tau)$ and increasing $a(\tau)$, ensuring good adhesion.

Equation (18) now integrates to give

$\tau=\frac{2}{|\lambda|^{2}}\left((X-1)|\lambda|+\log{\frac{1-|\lambda|X}{1-|\lambda|}}\right)\,,$
(19)

again determining $X$ implicitly as a function of $\tau$. Figure 6(b) shows
$X^{2}(\tau)$ as determined by (19) for several values of $|\lambda|$, both
less than and greater than 1. For $|\lambda|<1$, $X^{2}$ goes to zero at a
time $\tau^{*}(|\lambda|)$, indicating a finite-time singularity
[$h(\tau)=\infty$] at this time. If $|\lambda|=0$, the solution is

$X=(1-\tau)^{1/2}\quad\textnormal{so}\quad
h(\tau)=h_{0}\left(1-\tau\right)^{-1/4}\,,\quad
a(\tau)=a_{0}(1-\tau)^{1/8}\,,$ (20)

indicating the finite-time singularity at $\tau=1$ anticipated in the
introduction. These results hold only for so long as $h(\tau)\ll a(\tau)$, and
cease to apply as $\tau=t/|t_{0}|$ approaches the singularity time,
specifically when $1-\tau\sim(h_{0}/a_{0})^{8/3}$; at this stage, it is to be
expected that the two plates will rapidly separate. However, as will be shown
in §4 below, the evolution for $\tau>0$ is subject to a fingering instability.

More generally, if $F=F(t)$ is time-dependent, then (10a) still holds, giving

$\frac{1}{h^{4}}=\left(\frac{8\pi\,}{3\mu\,V^{2}}\right)\int_{0}^{t}F(t^{\prime})\textnormal{d}t^{\prime}\,+\frac{1}{h_{0}^{4}}.$
(21)

For prescribed $F(t)$, e.g. $F(t)\sim t^{\,\alpha}$, it is then a trivial
matter to obtain results analogous to (16); we shall not pursue these
possibilities here.

##  3 Effect of a trapped air bubble

Suppose now that a bubble of air is trapped so that the viscous liquid now
occupies the annular region $b(t)<r<a(t)$. The volume of the bubble is
$V_{b}=\pi b^{2}h$, and treating the air as incompressible (see §3.1 below),
this is constant, so that

$\frac{\textnormal{d}b}{\textnormal{d}t}=-\frac{b}{2h}\frac{\textnormal{d}h}{\textnormal{d}t}\,.$
(22)

The volume of the viscous liquid is now $V=\pi\left(a^{2}-b^{2}\right)h$, and
this is also constant. This leads to only minor changes in the above analysis.
Equation (4) obviously remains unchanged. Also the flux balance $2\pi r/h=-\pi
r^{2}\bar{u}(r,t)\textnormal{d}h/\textnormal{d}t$ still holds, since the air
as well as the viscous liquid is assumed incompressible. Hence (5) remains
valid also and the pressure is still given by (8) for $b<r<a$. The pressure
$p_{b}(t)$ in the air bubble is then given by

$p_{b}(t)=p(b,t)+\gamma/h=p_{a}-\frac{3\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\left(a^{2}-b^{2}\right)\,.$
(23)

The force $F=F_{p}$ is now given by

$F=\int_{b}^{\,a}\left(p(r,t)-p_{a}\right)\,2\pi
r\,\textnormal{d}r+(p_{b}-p_{a})\,\pi b^{2}\,.$ (24)

Using (8) and (23), this evaluates to

$F=-\frac{3\pi\mu}{2h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\left(a^{2}-b^{2}\right)\left(a^{2}+b^{2}\right)-\frac{\pi\gamma}{h^{2}}\left(a^{2}-b^{2}\right)=\frac{3\mu
V(V+2V_{b})}{8\pi}\frac{\textnormal{d}}{\textnormal{d}t}\left(\frac{1}{h^{4}}\right)-\frac{\gamma\,V}{h^{2}}\,,$
(25)

which reduces as expected to (10) when $V_{b}=0$.

Thus again, with $t_{0}$ now given by

$t_{0}=\frac{3\mu}{8\pi Fh_{0}^{4}}V(V+2V_{b})\,,$ (26)

we arrive at the following trivial modifications of (16) and (19), if surface
tension is neglected:

$\displaystyle h(\tau)$ $\displaystyle=$ $\displaystyle
h_{0}(1\\!+\\!\tau)^{-1/4}\,,\,a(\tau)=a_{0}(1\\!+\\!\tau)^{1/8}\,,\,b(\tau)=b_{0}(1\\!+\\!\tau)^{1/8}\,,\,\,\textnormal{if}\,\,F>0\,\textnormal{with}\,\tau=t/t_{0}\,,$
(27) $\displaystyle h(\tau)$ $\displaystyle=$ $\displaystyle
h_{0}(1\\!-\\!\tau)^{-1/4}\,,\,a(\tau)=a_{0}(1\\!-\\!\tau)^{1/8}\,,\,b(\tau)=b_{0}(1\\!-\\!\tau)^{1/8}\,,\,\,\textnormal{if}\,F<0\,\textnormal{with}\,\tau=t/|t_{0}|\,.$
(28)

Thus the presence of the bubble merely affects the time-scale of this
evolution, and the one-eighth power law still applies [now to both $a(\tau)$
and $b(\tau)$]. We shall find however that the stability of these evolving
states is seriously affected by the presence of the bubble.

###  3.1 Neglect of compressibility of air in the bubble

Note that, from (23), when $\textnormal{d}h/\textnormal{d}t<0$,
$p_{b}(t)>p_{a}$, i.e. the pressure in the bubble is greater than atmospheric
(this is essentially because $\bar{u}_{r}>0$, so that, according to (4),
$p(r,t)$ is monotonic decreasing in $r$); the air in the bubble is therefore
slightly compressed. This effect is however small in the experiments described
in §7 below; for, using (26) and (27), (23) becomes

$p_{b}(\tau)=p_{a}+\frac{2F\,h_{0}}{\left(V+V_{b}\right)(1+\tau)^{9/4}}\,.$
(29)

With the estimates $F\sim 10$ kg m s-2, $V+V_{b}\sim 10^{-5}$m3 and $h_{0}\sim
2\times 10^{-3}$m, we find that $2Fh_{0}/(V+V_{b})\sim 4\times 10^{3}$ kg
m-1s-2, small compared with atmospheric pressure $p_{a}\sim 10^{5}$ kg m-1s-2.
Thus at time $\tau$, $[p_{b}(\tau)-p_{a}]/p_{a}\sim 0.04(1+\tau)^{-9/4}$ and
Boyle’s Law implies a correspondingly small reduction in the volume of the
bubble due to compression. The compressibility of the air in the bubble is
therefore negligible in the conditions of our experiment, although other
circumstances can be imagined when the effect might become significant.

##  4 Stability of the basic state

In this section, we suppose $F>0$, and for simplicity we neglect the effect of
surface tension. We suppose first that the viscous liquid covers a simply-
connected domain $S\\!(t)$ in the $\\{x,y\\}$-plane, bounded by the curve
$C(t)\\!:\\!r=R(\theta,t)$, where $x=r\cos\theta,y=r\sin\theta$, and
$R(\theta,t)$ is single-valued and periodic in $\theta$, with period $2\pi$.
The area of $S\\!(t)$ is

$A(t)=\iint_{S}\textnormal{d}x\,\textnormal{d}y=\int_{0}^{2\pi}\left[\int_{0}^{R(\theta,t)}r\,\textnormal{d}r\right]\,\textnormal{d}\theta=\textstyle{\frac{1}{2}}\displaystyle\int_{0}^{2\pi}R^{2}(\theta,t)\,\textnormal{d}\theta\,,$
(30)

and the volume of the drop is $V=h(t)A(t)$. With origin $r=0$ at the centre-
of-mass of the drop,

$\iint_{S(t)}\\{x,y\\}\,\textnormal{d}x\,\textnormal{d}y=\textstyle{\frac{1}{3}}\displaystyle\int_{0}^{2\pi}R^{3}(\theta,t)\\{\cos\theta,\sin\theta\\}\,\textnormal{d}\theta=0\,.$
(31)

The pressure $p=p(r,\theta)$ still satisfies Reynolds’ equation (6) so that

$\nabla^{2}p=\frac{12\,\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\,,\quad
p=0\,\,\textnormal{on}\,\,C(t)\,.$ (32)

If $C(t)$ is known at any instant $t$, this Dirichlet problem may be solved
numerically; this provides the basis for a step-by-step determination of the
evolution of $C(t)$ from an initial condition $C(0)=C_{0}$, say.

###  4.1 Elliptic drop

Figure 8: Pressure contours $p=$ const for the particular initial condition
$a_{0}/b_{0}=3$, as given by (35); the magnitude of the pressure gradient at
the boundary is indicated by the length of the arrows.

The situation is well illustrated by the case of an elliptic drop, for which
$C(t)$ has the form

$\frac{x^{2}}{a^{2}(t)}+\frac{y^{2}}{b^{2}(t)}=1\,,\quad\textnormal{with}\quad
a(0)=a_{0},\,b(0)=b_{0},$ (33)

and where

$\pi\,a(t)b(t)h(t)=V=\textnormal{constant}=\pi\,a_{0}b_{0}h_{0}\,.$ (34)

We may assume $b_{0}\leq a_{0}$. The pressure is given by

$p(x,y,t)-p_{a}=\hat{p}(t)\left(\frac{x^{2}}{a^{2}(t)}+\frac{y^{2}}{b^{2}(t)}-1\right),$
(35)

so that, using (32),

$\nabla^{2}p=2\hat{p}(t)\left(\frac{1}{a^{2}(t)}+\frac{1}{b^{2}(t)}\right)=\frac{12\,\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\,.$
(36)

Hence

$\hat{p}(t)=\frac{6\,\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\frac{a^{2}b^{2}}{a^{2}+b^{2}}=\frac{6\,\mu\,V^{2}}{\pi^{2}h^{5}}\frac{\textnormal{d}h}{\textnormal{d}t}\frac{1}{a^{2}+b^{2}}\,.$
(37)

The downward force $F$ applied to the upper plate is now

$F=\iint_{S(t)}\hat{p}(t)\left(\frac{x^{2}}{a^{2}(t)}+\frac{y^{2}}{b^{2}(t)}-1\right)\,\textnormal{d}x\,\textnormal{d}y=-(1/2)\hat{p}(t)\,\pi
ab=-(V/2)\hat{p}(t)/h(t)\,.$ (38)

The pressure contours $p=$ const. are shown in figure 8. The pressure gradient
at the boundary, as shown by the arrows, is greater where these contours are
close together; the rate of expansion in the $y-$direction is therefore
greater than in the $x-$direction. We can calculate these rates as follows.

(a) (b)

Figure 9: Evolution of an elliptic drop, with initial condition
$a_{0}=1,b_{0}=0.5,\kappa=0.75$; (a) expanding drop when $F>0$, shown at times
$\tau=0,1,10,10^{2},10^{3}$ and $10^{4}$; the shape becomes gradually more
circular in form; (b) contracting drop when $F<0$; the minor axis of the
ellipse $b(\tau)$ decreases to zero at the singularity time
$\tau=0.71672356\dots$; the contracting boundary, shown at times
$\tau=0,0.6,0.7,0.7166,0.71672356$, is subject to smaller-scale instabilities.

The $x$-component of velocity at the point $(a,0)$ is

$\frac{\textnormal{d}a}{\textnormal{d}t}=\bar{u}(x,0)=-\frac{h^{2}}{12\mu}\left.\frac{\partial
p}{\partial x}\right|_{x=a,\,y=0}=-\frac{h^{2}}{6\mu
a}\hat{p}(t)=-\frac{1}{a}\frac{V^{2}}{\pi^{2}h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\frac{1}{(a^{2}+b^{2})}$
(39)

and hence

$\frac{\textnormal{d}a^{2}}{\textnormal{d}t}=\frac{V^{2}}{\pi^{2}}\frac{\textnormal{d}}{\textnormal{d}t}\left(\frac{1}{h^{2}}\right)\frac{1}{\left(a^{2}+b^{2}\right)}\,.$
(40)

Similarly,

$\frac{\textnormal{d}b^{2}}{\textnormal{d}t}=\frac{V^{2}}{\pi^{2}}\frac{\textnormal{d}}{\textnormal{d}t}\left(\frac{1}{h^{2}}\right)\frac{1}{\left(a^{2}+b^{2}\right)}$
(41)

and it follows that $a^{2}(t)-b^{2}(t)=$ constant $=a_{0}^{2}-b_{0}^{2}$. We
then have from (40)

$\frac{\textnormal{d}}{\textnormal{d}t}\left[a^{4}-\left(a_{0}^{2}-b_{0}^{2}\right)a^{2}\right]=\frac{V^{2}}{\pi^{2}}\frac{\textnormal{d}}{\textnormal{d}t}\left(\frac{1}{h^{2}}\right),$
(42)

so that, with the initial condition $a(0)=a_{0}$, $a^{2}$ satisfies the
quadratic equation

$a^{4}-(a_{0}^{2}-b_{0}^{2})a^{2}=V^{2}/\pi^{2}h^{2}\,,$ (43)

with relevant root

$a^{2}(t)=\textstyle{\frac{1}{2}}\displaystyle\left(a_{0}^{2}-b_{0}^{2}\right)+\textstyle{\frac{1}{2}}\displaystyle\left[\left(a_{0}^{2}-b_{0}^{2}\right)^{2}+4V^{2}/\pi^{2}h^{2}\right]^{1/2}\,.$
(44)

The corresponding expression for $b^{2}(t)$ is

$b^{2}(t)=a^{2}(t)-\left(a_{0}^{2}-b_{0}^{2}\right)=-\textstyle{\frac{1}{2}}\displaystyle\left(a_{0}^{2}-b_{0}^{2}\right)+\textstyle{\frac{1}{2}}\displaystyle\left[\left(a_{0}^{2}-b_{0}^{2}\right)^{2}+4V^{2}/\pi^{2}h^{2}\right]^{1/2}\,,$
(45)

so that

$a^{2}(t)+b^{2}(t)=\left[\left(a_{0}^{2}-b_{0}^{2}\right)^{2}+4V^{2}/\pi^{2}h^{2}\right]^{1/2}\,.$
(46)

From (37) we then have

$\hat{p}(t)=-\frac{6\,\mu\,V^{2}}{\pi
h^{4}}\frac{\textnormal{d}h}{\textnormal{d}t}\left[\left(a_{0}^{2}-b_{0}^{2}\right)^{2}\pi^{2}h^{2}+4V^{2}\right]^{-1/2}\,.$
(47)

and from (48),

$F=\frac{3\,\mu\,V^{2}}{8\pi}\frac{\textnormal{d}}{\textnormal{d}t}\left(\frac{1}{h^{4}}\right)\left[1+\left\\{\left(a_{0}^{2}-b_{0}^{2}\right)\pi
h/2V\right\\}^{2}\right]^{-1/2}\,,$ (48)

in agreement with (10) when $a_{0}=b_{0}$.

In dimensionless form, with $\tau=t/t_{0}$ [where $t_{0}$ is still defined by
(11)] and $\hat{h}(\tau)=h(t)/h_{0}$, (48) becomes

$\frac{\textnormal{d}}{\textnormal{d}\tau}\left(\frac{1}{\hat{h}^{4}}\right)=\left[1+\kappa^{2}\hat{h}^{2}\right]^{1/2}\,,$
(49)

where

$\kappa=\left(a_{0}^{2}-b_{0}^{2}\right)/2a_{0}b_{0}\,.$ (50)

Here, $\kappa\geq 0$ is a measure of the initial eccentricity of the ellipse.
With this notation, (44) and (45) become

$a^{2}(\tau)/a_{0}b_{0}=\kappa+\left[\kappa^{2}+\hat{h}^{-2}(\tau)\right]^{1/2}\,,\quad
b^{2}(\tau)/a_{0}b_{0}=-\kappa+\left[\kappa^{2}+\hat{h}^{-2}(\tau)\right]^{1/2}\,,$
(51)

The evolution is shown in figure 9 for initial conditions $a_{0}=1,b_{0}=0.5$
(so $\kappa=0.75$). In (a), $F>0$, and the drop expands, becoming gradually
more circular in form. In (b), $F<0$, and the drop contracts with increasing
ellipticity [as $\hat{h}(\tau)$ increases], becoming singular as
$\tau\rightarrow\tau_{c}=0.71672\dots\,$, when $b(\tau)\rightarrow 0$ and
$a(t)\rightarrow\surd{3}/2=0.8660\dots\,$. However, this behaviour persists
only for so long as $b(\tau)\gtrsim h_{0}\,\hat{h}(\tau)$, and it is moreover
subject to a ‘fingering instability’. Surface tension must also become
important where the curvature of the boundary becomes large. These effects
will be considered in §5 below.

###  4.2 Linearised stability analysis

We here consider a circular drop with perturbed boundary

$R(\theta,0)=a_{0}+\epsilon\,\alpha_{0}\cos n\theta$ (52)

at the initial instant, where $n$ is a positive integer. When $n=1$, the
circle is merely displaced in the $x$-direction without distortion; we may
therefore suppose that $n\geq 2$, so that (31) is automatically satisfied.
(When $n=2$, the circle is perturbed to an ellipse.) With the choice (52),

$A(0)=\pi
a_{0}^{2}+\textstyle{\frac{1}{2}}\displaystyle\pi\,\epsilon^{2}\alpha_{0}^{2}.$
(53)

When $0<\epsilon\ll 1$, the problem may be linearised. Then for $t>0$,

$R(\theta,t)=a(t)+\epsilon\,\alpha(t)\cos n\theta\,,$ (54)

where $a(t)$ is given by (16) and $A(t)=\pi
a^{2}(t)+\textnormal{O}(\epsilon^{2})$. Thus, let

$p(r,\theta,t)=p_{0}(r,t)+\epsilon\,p_{1}(r,t)\cos n\theta\,,$ (55)

where $p_{0}(r,t)$ is given by (8). At this level of approximation, $h(t)$ is
still given by (16), and since
$\nabla^{2}p_{0}=12\,\mu\,h^{-3}\textnormal{d}h/\textnormal{d}t$, (6) gives

$\frac{1}{r}\frac{\partial}{\partial r}r\frac{\partial p_{1}}{\partial
r}-\frac{n^{2}}{r^{2}}p_{1}=0\,.$ (56)

The solution, finite at $r=0$, is

$p_{1}(r,t)=k(t)r^{n}\,,$ (57)

so that now

$p(r,\theta,t)=p_{0}(r,t)+\epsilon\,k(t)r^{n}\cos n\theta\,.$ (58)

The pressure boundary condition is now $p(r,\theta,t)-p_{a}=0$ on
$r=a(t)+\epsilon\,\alpha\cos n\theta$, and, noting that

$p_{0}(a(t)+\epsilon\,\alpha(t)\cos
n\theta,t)-p_{a}=p_{0}(a)-p_{a}+\epsilon\,\alpha(t)\cos
n\theta\left.\frac{\partial p_{0}}{\partial
r}\right|_{r=a}=\epsilon\,a(t)\,\alpha(t)\cos
n\theta\frac{6\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}$ (59)

and that

$\epsilon\,p_{1}(a(t)+\epsilon\,\alpha(t)\cos
n\theta,t)=\epsilon\,k(t)a^{n}+\textnormal{O}(\epsilon^{2})\,,$ (60)

this boundary condition (at order $\epsilon$) determines $k(t)$, leading
finally to

$p(r,\theta,t)-p_{a}=\left[(r^{2}-a^{2})-2\epsilon\,a\,(r/a)^{n}\,\alpha(t)\cos
n\theta\right]\frac{3\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}+\textnormal{O}(\epsilon^{2})\,,$
(61)

and so

$\frac{\partial p}{\partial r}=\left[r-\epsilon\,n\,\alpha(t)(r/a)^{n-1}\cos
n\theta\right]\frac{6\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}+\textnormal{O}(\epsilon^{2})\,.$
(62)

On $r=R(\theta,t)=a+\epsilon\,\alpha\cos n\theta$, this therefore gives

$\frac{\partial p}{\partial r}=[a-\epsilon\,(n-1)\,\alpha\cos
n\theta]\frac{6\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}+\textnormal{O}(\epsilon^{2})\,.$
(63)

Now the curve $C(t)$ moves according to the equation

$\frac{\partial R}{\partial
t}=\frac{\textnormal{d}a}{\textnormal{d}t}+\epsilon\,\frac{\textnormal{d}\alpha}{\textnormal{d}t}\cos
n\theta=\bar{u}_{r}(R,t)+\textnormal{O}(\epsilon^{2})\,,$ (64)

and $\bar{u}_{r}$ still satisfies (4), i.e.

$\bar{u}_{r}=-\frac{h^{2}}{12\mu}\,\left.\frac{\partial p}{\partial
r}\right|_{r=a(t)+\epsilon\,\alpha(t)\cos n\theta}\,.$ (65)

At leading order, (64) merely confirms that $a^{2}(t)\,h(t)=$ const. At order
$\epsilon$, with $h(t)$ given by (16), we obtain the equation for $\alpha(t)$,

$\frac{\textnormal{d}\alpha}{\textnormal{d}t}=\frac{(n-1)}{2h}\frac{\textnormal{d}h}{\textnormal{d}t}\,\alpha=-\frac{(n-1)\,\alpha}{8(t+t_{0})}\,,$
(66)

which integrates to give

$\alpha(t)=\alpha_{0}(1+t/t_{0})^{-(n-1)/8}.$ (67)

Since $n\\!>\\!1$, this ensures that, when $F\\!>\\!0$ and so $t_{0}\\!>\\!0$,
the perturbation decays to zero as $t\rightarrow\infty$; and in fact the
higher the value of $n$, the more rapid is this decay.

(a) (b)

Figure 10: (a) Sketch indicating the stabilising mechanism for a squeezed
drop; the drop expands more rapidly where the isobars $p=$ const. are closer
together, thus tending to restore a circular shape; (b) sketch indicating the
destabilising mechanism at the boundary of an air bubble, in which the
pressure $p_{b}$ is greater than atmospheric pressure $p_{a}$; the expansion
is again more rapid where the isobars are closer together, leading now to
growth of the perturbation amplitude.

The stability can be best understood from consideration of the sketch of
figure 10(a) which shows the particular situation when $n=5$ (a five-fold
hypotrochoid). Just as for the case of an elliptic drop (figure 8), the
restoring force is greater where the isobars are closer, thus tending to
restore a circular shape. It is here that our conclusion differs from that of
Savina & Nepomnyashchy (2015), as mentioned in the introduction.

##  5 Fingering instability when $F\\!<0$

Figure 11: Plot of $\hat{\alpha}(\tau)=\alpha(\tau)/\alpha_{0}$ as a function
of $\tau=t/|t_{0}|$, for $\sigma\,|t_{0}|=0.01$ and $n=2,8,12,15$.

When $F\\!<0$, the basic state is given by (19), so that, for $0<t\ll|t_{0}|$,

$\frac{1}{h}\frac{\textnormal{d}h}{\textnormal{d}t}=\frac{1}{4(|t_{0}|-t)}\,,$
(68)

and (66) becomes

$\frac{\textnormal{d}\alpha}{\textnormal{d}t}=\frac{(n-1)\,\alpha}{8(|t_{0}|-t)}\,,$
(69)

which integrates to give

$\alpha(t)=\alpha_{0}(1-t/|t_{0}|)^{-(n-1)/8}\,.$ (70)

This shows the beginning of a ‘fingering instability’ of a type first
identified by Saffman & Taylor (1958), resulting here from the radial
contraction of the region occupied by the viscous liquid and the resulting
suction of air into this region. The instability is stronger as $n$ increases,
apparently without limit as $n\rightarrow\infty$. In reality, surface tension
must play a role in determining the most unstable mode, an effect that may be
simply treated as follows.

Let $\gamma$ be the surface tension acting at the boundary
$r=R(\theta,t)=a+\epsilon\,\alpha\cos n\theta$. As we have seen in §2, this
leads to a jump in pressure $\gamma/h(t)$ across the air/liquid interface;
this jump, independent of $\theta$, is irrelevant for the stability of the
circular shape and can be ignored here. There is also however a contribution
to the pressure jump from the variation of curvature around the drop boundary.
We assume that $\gamma$ is relatively weak, so that surface tension is
significant only for modes $\sim\\!\cos n\theta$ for which $n\gg 1$. Then the
pressure condition $p(r,t)-p_{a}=0$ on $r=R$ should be replaced by

$p(r,t)-p_{a}\approx-\gamma\frac{1}{a^{2}}\frac{\partial^{2}R}{\partial\theta^{2}}=\gamma
n^{2}\frac{\alpha}{a^{2}}\cos n\theta\quad\textnormal{on}\quad
r=a+\epsilon\,\alpha\cos n\theta\,.$ (71)

This modified boundary condition leads to a corresponding modification in
(69), which now becomes, after some simplification,

$\frac{\textnormal{d}\alpha}{\textnormal{d}t}=\frac{(n-1)\,\alpha}{8|t_{0}|(1-t/|t_{0}|)}-\frac{\sigma\,n^{3}\alpha}{8(1-t/|t_{0}|)^{7/8}}\quad\textnormal{where}\quad\sigma=\left(\frac{2\gamma}{3\mu}\right)\frac{h_{0}^{2}}{a_{0}^{3}}\,.$
(72)

With $\tau=t/|t_{0}|$, this integrates to give

$\alpha(\tau)=\alpha_{0}(1-\tau)^{-(n-1)/8}\exp{\left[-n^{3}\sigma\,|t_{0}|\left\\{1-(1-\tau)^{1/8}\right\\}\right]}\,.$
(73)

By way of example, the function $\alpha(\tau)/\alpha_{0}$ is shown in figure
11, for the parameter values $\sigma\,|t_{0}|=0.01$ and $n=2,8,12$ and $15$,
and for the interval $0<\tau<0.99$. The modes become increasingly unstable as
$n$ increases from 2 to 8, but by $n=12$, the stabilising effect of surface
tension is evident, and the modes are rapidly damped for $n\gtrsim 15$. More
generally, surface tension has a strong stabilising effect for modes for which
$n\gtrsim(\sigma\,|t_{0}|)^{-1/3}$.

##  6 Destabilising effect of trapped air bubble

We here follow the procedure described in §4.2 above, to examine the stability
of the evolving state with a trapped air bubble, as described in §3. For
simplicity in this section, we ignore surface tension and assume first that
$F>0$ so that $\textnormal{d}h/\textnormal{d}t<0$. In the undisturbed state,
the viscous liquid is contained in the region $a(t)>r>b(t)$. We suppose that
the liquid boundaries are perturbed to

$r=R_{1}(\theta,t)=a(t)+\epsilon\,\alpha(t)\cos{n\theta}\,,\quad
r=R_{2}(\theta,t)=b(t)+\epsilon\,\beta(t)\cos{n\theta}\,,$ (74)

and linearise in $\epsilon$. The pressure field in the region
$R_{1}(\theta,t)>r>R_{2}(\theta,t)$ still takes the form

$p(r,\theta,t)=p_{a}+p_{0}(r,t)+\epsilon p_{1}(r,\theta,t)\,,$ (75)

where $\nabla^{2}p_{1}=0$, but now

$p_{1}(r,\theta,t)=-\frac{6\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\,\left[\hat{k}(t)r^{n}+\hat{m}(t)r^{-n}\right]\cos{n\theta}\,,$
(76)

where $\hat{k}(t)$ and $\hat{m}(t)$ are to be determined. As before, the
condition $p(r,\theta,t)=p_{a}$ on $r=R_{1}(\theta,t)$ leads, at order
$\epsilon$, to

$\hat{k}(t)\,a^{n}(t)+\hat{m}(t)\,a^{-n}(t)=\alpha(t)\,a(t)$ (77)

and similarly the condition $p(r,\theta,t)-p_{a}=p_{b}$ on $r=R_{2}(\theta,t)$
leads, at order $\epsilon$, to

$\hat{k}(t)\,b^{n}(t)+\hat{m}(t)\,b^{-n}(t)=\beta(t)\,b(t)\,.$ (78)

Solving (77) and (78) for $\hat{k}(t)$ and $\hat{m}(t)$, we obtain

$\hat{k}(t)=\frac{\alpha\,a^{n+1}-\beta\,b^{n+1}}{a^{2n}-b^{2n}}\qquad\textnormal{and}\quad\hat{m}(t)=\frac{a^{n}b^{n}\left(\beta\,b\,a^{n}-\alpha\,a\,b^{n}\right)}{a^{2n}-b^{2n}}\,.$
(79)

We now have

$\frac{\partial p}{\partial
r}=\frac{6\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\,r+\epsilon\left[n\,\hat{k}\,r^{n-1}-n\,\hat{m}\,r^{-n-1}\right]\cos{n\theta}\,,$
(80)

so that, on $r=R_{2}=b+\epsilon\,\beta\cos{n\theta}$,

$\left.\frac{\partial p}{\partial
r}\right|_{r=R_{2}}=\frac{6\mu}{h^{3}}\frac{\textnormal{d}h}{\textnormal{d}t}\,(b+\epsilon\,\beta\cos{n\theta})+\epsilon\left[n\,\hat{k}\,b^{n-1}-n\,\hat{m}\,b^{-n-1}\right]\cos{n\theta}+\textnormal{O}(\epsilon^{2})\,.$
(81)

Hence, recalling (4), we have

$\left.\bar{u}_{r}\right|_{r=R_{2}}=\left.\bar{u}_{r}\right|_{0}-\frac{1}{2h}\frac{\textnormal{d}h}{\textnormal{d}t}\,\left[\beta+n\left(\hat{k}\,b^{n-1}-\hat{m}\,b^{-n-1}\,\right)\right]\epsilon\,\cos{n\theta}\,,$
(82)

where $\left.\bar{u}_{r}\right|_{0}$ is the velocity in the unperturbed state.
At order $\epsilon$, $\partial R_{2}/\partial
t=\left.\bar{u}_{r}\right|_{r=R_{2}}$ now leads to

$\frac{\textnormal{d}\beta}{\textnormal{d}t}=-\frac{1}{2h}\frac{\textnormal{d}h}{\textnormal{d}t}\,\left\\{\beta-n\left[\frac{2a^{n+1}b^{n-1}\alpha-\left(a^{2n}+b^{2n}\right)\beta}{a^{2n}-b^{2n}}\right]\right\\}\,.$
(83)

Similarly, $\partial R_{1}/\partial t=\left.\bar{u}_{r}\right|_{r=R_{1}}$
leads at order $\epsilon$ to

$\frac{\textnormal{d}\alpha}{\textnormal{d}t}=-\frac{1}{2h}\frac{\textnormal{d}h}{\textnormal{d}t}\,\left\\{\alpha-n\left[\frac{\left(a^{2n}+b^{2n}\right)\alpha-2a^{n-1}b^{n+1}\beta}{a^{2n}-b^{2n}}\right]\right\\}\,.$
(84)

With dimensionless time $\tau=t/t_{0}$, noting from (27) that

$-\frac{1}{2h}\frac{\textnormal{d}h}{\textnormal{d}t}=\frac{1}{8(1+\tau)}\quad\textnormal{and}\quad\frac{b(\tau)}{a(\tau)}=\textnormal{const.}=\eta,\,\textnormal{say,}\quad\textnormal{where}\,\,0<\eta<1\,,$
(85)

these equations may be written in the form

$\frac{\textnormal{d}\alpha}{\textnormal{d}\tau}=\frac{1}{8(1+\tau)}(A\,\alpha+B\,\beta)\,,\quad\frac{\textnormal{d}\beta}{\textnormal{d}\tau}=\frac{1}{8(1+\tau)}(C\,\alpha+D\,\beta)\,,$
(86)

where

$A=1-n\left(\frac{1+\eta^{2n}}{1-\eta^{2n}}\right)\,,\quad
B=\frac{2n\,\eta^{n+1}}{1-\eta^{2n}}\,,\quad
C=-\frac{2n\,\eta^{n-1}}{1-\eta^{2n}}\,,\quad
D=1+n\left(\frac{1+\eta^{2n}}{1-\eta^{2n}}\right)\,.$ (87)

Note that

$A+D=2,\quad\textnormal{and}\quad AD-BC=1-n^{2}\,.$ (88)

The equations (86) admit solutions of the form
$(\alpha,\beta)=\left(\hat{\alpha},\hat{\beta}\right)\,(1+\tau)^{\sigma}$,
provided

$8\sigma\,\hat{\alpha}=A\,\hat{\alpha}+B\,\hat{\beta}\,,\quad
8\sigma\,\hat{\beta}=C\,\hat{\alpha}+D\,\hat{\beta}\,.$ (89)

The determinant condition for a non-trivial solution reduces to

$8\sigma=\textstyle{\frac{1}{2}}\displaystyle(A+D)\pm\sqrt{1-(AD-BC)}=1\pm
n\,,$ (90)

using the results (88). The modes for which $\sigma=\sigma_{1}=(1+n)/8$ and
$\sigma=\sigma_{2}=(1-n)/8$ are respectively unstable and stable, and the
corresponding ratios $\hat{\alpha}/\hat{\beta}$ for these modes simplify to

$\left[\hat{\alpha}/\hat{\beta}\right]_{1}=\eta^{n+1}\quad\textnormal{and}\quad\left[\hat{\alpha}/\hat{\beta}\right]_{2}=\eta^{1-n}.$
(91)

For example, if $\eta=\textstyle{\frac{1}{2}}\displaystyle$ and $n=5$, then
$\left[\hat{\alpha}/\hat{\beta}\right]_{1}=1/64$ and
$\left[\hat{\alpha}/\hat{\beta}\right]_{2}=16$. Thus, the unstable mode has
much larger amplitude at the bubble boundary $r=b$, whereas the stable mode
has much larger amplitude at the outer boundary $r=a$; we may think of the
instability as being ‘driven’ from the bubble boundary. These properties
become stronger as $\eta$ decreases and/or $n$ increases.

Figure 12: Solution of equations (86) for $\alpha(\tau)$ and $\beta(\tau)$
with initial conditions $\alpha(0)=\beta(0)=1$ and with $\eta=1/2,\,n=5$, as
in the sketch of figure 10(b); both $\alpha(\tau)$ and $\beta(\tau)$
ultimately increase at the same rate, the ratio $\beta(\tau)/\alpha(\tau)$
tending to a constant for large $\tau$. This behaviour is typical.

The general solution corresponding to arbitrary initial conditions
$\alpha(0)=\alpha_{0},\beta(0)=\beta_{0}$ is of course a superposition of
these modal solutions, but the unstable mode rapidly dominates. The mechanism
of this instability is indicated in the sketch of figure 10(b): the air bubble
is expanding in the mean, but the expansion is more rapid where the pressure
contours (isobars) are closer together; this obviously leads to increase of
the perturbation amplitude.

Figure 12 shows the behaviour when $\alpha(0)=\beta(0)=1,\eta=b/a=1/2$ and
$n=5$. Here, the unstable mode, for which $\beta(\tau)/\alpha(\tau)=64$,
dominates for $\tau\gtrsim 1.4$. The growth rate of this type of instability
increases with $n$, and is controlled for large $n$ by surface tension, just
as in §5; we need not labour the details here.

###  6.1 The situation when $F<0$

When $F<0$, with the basic state given by (28) there are still two
perturbation modes, one stable and one unstable, but the unstable mode is now
‘driven’ by the fingering instability at the outer boundary $r=a(t)$, where
the amplitude is much greater.

##  7 Simple experimental demonstrations

###  7.1 Single expanding drop

As already mentioned, our ‘home experiments’ were carried out under
effectively lockdown conditions.The first experiment was designed to verify
the one-eighth power law (16). For fluid we used Lyle’s Black
Treacle111http://fiches.ranson.be/00006253nl-be.pdf (density 1.41 g/ml,
viscosity $\mu\approx 64.3\,\textnormal{kg}\textnormal{
s}^{-1}\textnormal{m}^{-1}$ at 22oC). A drop of volume $V=5$ ml was placed on
a horizontal glass plate and a second glass plate of mass $M=1.127$ kg was
placed on this and allowed to settle, so that the downward force was
$F=Mg=11.04\,{\textnormal{kg}\,\textnormal{m}\,\textnormal{s}^{-2}}$; thus
$\mu/F=5.82\,\textnormal{m}^{-2}\textnormal{s}$. The plates were illuminated
from below, and a camera was set to take shots from above at 10-second
intervals. the initial time $t=0$ was set when the drop radius was
$a_{0}=30.72$ mm; thus $h_{0}=V/\pi a_{0}^{2}=1.69$ mm, and
$h_{0}/a_{0}=0.055$. From (11), the time-scale $t_{0}$ for this run was thus

$t_{0}=\frac{3\pi^{3}\,(5.82)(0.03072)^{8}}{8(5\times
10^{-6})^{2}}=2.15\,\textnormal{s}\,.$ (92)

Figure 1 shows photos of the expanding drop at three times $\tau=t/t_{0}=$ 60,
172, and 321. The background grid allowed measurement of the diameter in each
case. Figure 7(a) shows the growth of the drop radius for $0\leq\tau\leq 321$.
The measurements agree well with the result (16) as shown by the red curve;
and figure 7(b) shows that the one-eighth power law is remarkably well
satisfied over this full time range.

###  7.2 Contraction of drop by leverage

Figure 2 shows the effect of gently levering the plates open at the lower
right-hand corner some time after the stage of figure 1(c) was reached. The
fingering instability immediately develops and the fingers grow as the
levering is gradually increased. The duration of course depends on the speed
at which the angle between the plates is increased, and as the levering was
applied by hand, exact repetition was not possible. This particular sequence
was taken over a period of approximately one minute, and is quite typical. The
characteristic tip-splitting and side-branching of fingers is evident, and
cavitation bubbles can be detected in the very low pressure region in panel
(e). The reason for this cavitation may be understood with reference to the
‘hinged-plate’ problem briefly considered in §2 of Moffatt (1964); if two
plates $\theta=\pm\alpha(t)$ are hinged at $r=0$ with viscous fluid in the gap
$-\alpha(t)<\theta<\alpha(t)$, and if $\alpha(t)$ is slowly increased from a
very small value, then the pressure $p(r,t)$ behaves like
$\mu\,\dot{\alpha}(t)\log r$ which becomes negative as $r$ decreases towards
zero, so that cavitation is inevitable in this region. This is of course a
two-dimensional idealisation, but the same physical mechanism is presumably
responsible for the cavitation bubbles appearing in the present experiment.

We should note further that, as a consequence of the no-slip condition, a thin
residual film is left on both plates in the regions invaded by the fingers.
The thickness $s$ of these residual films is presumably given by $s/h\sim\mu
U/\gamma$, the result obtained by Bretherton (1961) for a bubble rising in a
capillary tube containing viscous liquid. In our present situation, with
$\mu\sim 70\,\textnormal{kg}\textnormal{ s}^{-1}\textnormal{m}^{-1}$,
$\gamma\sim 0.07$ kg s-2 and the observed creeping velocity $U$ in the range
$10^{-4}$–$10^{-5}$m s-1, this gives $s/h\sim 10^{-1}$–$10^{-2}$. When $h\sim
0.1$mm, the residual film thickness $s$ is therefore in the range 1–10 $\mu$m.

###  7.3 Evolution of an annular drop

Figure 3 shows an evolution sequence when the initial shape of the drop [panel
(a)] is as near to a uniform circular annulus as could be achieved by hand.
The curvature of the inner boundary of the annulus was inevitably not quite
uniform, and where it was maximal a perturbation immediately developed, as
anticipated by the stability analysis of §6. This perturbation evolves into a
secondary bubble [panel (b)] which is ejected, leaving the primary bubble
[panel (c)] which deforms slowly to a distinctly non-circular shape
[panel(d)]; the outer boundary was only weakly perturbed from a circle, The
secondary bubble detaches from the primary bubble and erupts through the outer
boundary [panels (d,e)] while a new (tertiary) bubble grows very slowly from
the residual cusp on the primary bubble. After about 45 minutes, the situation
[panel(f)] appears to be quasi-static, and very small wrinkles appear on the
outer boundary, probably associated with slight roughness of the two plates.
After a further hour, leverage was introduced at the lower right-hand corner,
and fingering rapidly developed [panels (g-i)] over about one minute. A
growing finger punctures the tertiary bubble [panel (g)] at the point marked
by the arrow, providing a connected path from the primary bubble to the
exterior. Fingers also erupt from the bubble [panel (j)] and encounter the
‘invading’ fingers, forming a defensive ridge. The pattern continues to evolve
round the outer boundary [panel (k)], and cavitation bubbles may be detected
in the ‘north-western’ sector ahead of the advancing fronts. Panel [l] shows
the ultimate fingering pattern just before rupture of the film when the plates
separate completely; the ridge, which remains a prominent feature of this
pattern, is breached only at the location of the puncture point.

###  7.4 Merging of four drops

We then carried out an experiment starting with four drops placed as near to
the corners of a square as could be achieved, as shown in figure 4 [panel
(a)]. As the upper plate descends, these drops expand and ultimately merge,
trapping an air bubble [panel (e)]. The initial drop contact [panel (c)]
involves an instantaneous geometrical singularity at a time $\tau=\tau_{c}$,
which is subsequently (i.e. for $\tau>\tau_{c}$) resolved by surface tension.
[We should distinguish such a singularity from the more familiar finite-time
singularities, as treated for example by Eggers & Fontelos (2015), for which a
singularity develops as $\tau\rightarrow\tau_{c}$ and is resolved by surface
tension or some other mechanism just before the singularity occurs, i.e. for
$\tau<\tau_{c}$.] After the bubble is trapped, it shows no tendency to become
circular; on the contrary, where the curvature on the air bubble is maximal, a
secondary bubble begins to pinch off [panel (g)], again leaving the very thin
residual films of treacle on both plates.

The formation of the secondary bubble is consistent with the qualitative
description of the bubble instability in §6 above: the pressure contours are
compressed where the bubble curvature is large, leading to amplification of
the perturbation in this region. The secondary bubble slowly propagates
towards the outer boundary, separates from the primary bubble (a singular
process analogous to the breaking of a viscous thread), and ultimately erupts
through this boundary [panel (i)] – again a singular process. Up to this
point, the total volume of trapped air is constant, but this falls very
rapidly at the moment of eruption. At the same time, a tertiary bubble begins
to form under the continuing influence of the pressure gradient; this suggests
an iterative process, although much slower at each stage.

###  7.5 Contraction by levering

Figure 5 shows a continuation of the sequence of figure 4, showing first
[panels (a-c)] the eruption of the secondary bubble and the very slow
separation of the tertiary bubble from the primary bubble which is now much
reduced in size. Note that the outer boundary remains relatively smooth but
does not become circular; this is because the internal instability has a
persistent weak disturbing effect on the outer boundary [see the comment
following (91)]. In panel (d), levering has just commenced at the lower right-
hand corner of the plates, and the first signs of the fingering instability
are visible. In panel (e), fingers erupting from the outer boundary of the
drop are evident, and tip-splitting and side-branching of the invading fingers
are well developed; in panel (f) a ridge develops between the invading and
defending cohorts of fingers. In panel (g) the defending cohort, now in
retreat, is punctured by an invading finger; the pressure is now atmospheric
($p_{a}$) throughout the fingers and the bubbles, but less than atmospheric in
the remaining bulk of treacle; in panel (h) the invading front of fingers
advances round the boundary and cavitation bubbles appear ahead of the front.
Finally, panel (i) shows that fingering has again extended throughout the
drop, to the extent that the treacle is almost entirely contained in the
narrow ‘tree branches and twigs’ that separate the fingers. The ultimate
pattern retains an imprint of the original trapped bubbles, in the pronounced
pentagonal ‘ridge’, still breached at only one point in the south-east sector.

##  8 Conclusions

The main conclusions of this investigation may be summarised as follows.
(i) A single drop of viscous fluid subject to squeezing under a constant force
$F\\!>\\!0$ in a Hele-Shaw cell tends to adopt a circular form, with radius
$a(t)$ increasing like $t^{1/8}$ (the ‘one-eighth power law’); this is
confirmed experimentally.
(ii) An initially elliptic drop subject to similar squeezing tends to circular
form; more generally, any perturbation of an expanding drop of circular form
decreases algebraically in time, i.e. the circular form is stable.
(iii) If the same drop is subject to contraction $(F\\!<0)$, then it is
unstable to a fingering disturbance whose scale is determined by surface
tension; this fingering establishes a lattice that ultimately extends
throughout the fluid domain, before rupture of the film and complete
separation of the plates.
(iv) Sufficiently large surface tension can lead to growth of $a(t)$ even when
$F\\!<0$, an important aspect of adhesion dynamics.
(v) An annular drop subject to squeezing can also follow a one-eighth power
law expansion of both inner and outer boundaries; but this expansion is
unstable, the instability growing predominantly from the inner boundary which
shows no tendency to remain circular.
(vi) The instability manifests itself through ejection of a secondary bubble
from the primary bubble leaving a residual layer on both plates whose
thickness is determined by surface tension; this secondary bubble erupts
through the outer boundary of the annular drop, and a tertiary bubble emerges
from the residual cusp on the primary bubble; this process is extremely slow.
(vii) When levering is applied to separate the plates, the fingering
instability occurs at the outer boundary of the annulus, and rapidly spreads
towards the inner boundary where much weaker fingers emerge; the bubble is
punctured at a single point by an invading finger, leading to equalisation of
the pressure in the bubble with the external atmosphere. It is like the
puncturing of an inflated balloon.
(viii) As fingering continues to spread round the outer boundary, cavitation
bubbles appear in the sector opposite the point of leverage, where the
pressure would otherwise fall below the vapour pressure of the treacle;
fingering ultimately extends over the whole fluid domain just before the final
sudden rupture of the film. A ridge is formed where the invading cohort of
fingers from the outer boundary impacts the defending cohort of fingers that
emerge from the bubble.
(ix) A similar sequence of events occur when the air bubble is trapped by the
expansion of four initially disjoint drops placed at the corners of a square;
in this case the cusp singularities that occur when the drops make contact are
resolved by surface tension as the drops continue to expand. Under subsequent
leverage, the processes of fingering, puncturing, ridge formation and
cavitation occur in similar manner in this case also.

## References

  * Ball & Huppert (2019) Ball, T.V. & Huppert, H.E. 2019 Similarity solutions and viscous gravity current adjustment times. J. Fluid Mech. 874, 285-298.
  * Bikerman (1947) Bikerman, J. J.1947 The fundamentals of tackiness and adhesion. J. Colloid Sci.  2, 163–175.
  * Bretherton (1961) Bretherton, F. P. 1961 The motion of long bubbles in tubes. J.Fluid Mech.  10, 166–188.
  * Davis (1974) Davis, P. J. 1974 The Schwarz Function and its Applications, The Carus Mathematical Monographs, No. 17 (Math. Assoc. of America) 228 pp.
  * Eggers & Fontelos (2015) Eggers, J. & Fontelos, M.A. 2015 Singulariities: Formation, Structure, and Propagation. Cambridge Texts in Applied Mathematics, Cambridge University Press, 453 pp.
  * Gaudet, McKinley & Stone (1996) Gaudet, S., McKinley, G. H. & Stone, H.A. 1996 Extensional deformation of Newtonian liquid bridges. Phys. Fluids, 8, 2568–2579.
  * Gay (2002) Gay, C. 2002 Stickiness – some fundamentals of adhesion. Integr. Compar. Biol., 42, 1123–1126.
  * Gustafsson & Vasiliev (2006) Gustafsson, B. & Vasiliev, A. 2006 Conformal and Potential Analysis in Hele-Shaw Cells, Birkhäuser Verlag, Basel-Boston-Berlin, 183 pp.
  * Howison (2000) Howison, S.D. 2000 A note on the two-phase Hele-Shaw problem. J.Fluid Mech.  409, 243–249.
  * Huppert (1982) Huppert, H.E. 1982 The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface. J.Fluid Mech.  121, 43–58.
  * Huppert et al. (1982) Huppert, H.E., Shepherd, J.B., Sigurdsson, H & Sparks, S.J. 1982 On lava dome growth, with application to the 1979 lava extrusion of the Soufrière of St Vincent. J. Volcan. Geoth. Res., 14, 199–222.
  * Kelly & Hinch (1997) Kelly, E.D. & Hinch, E.J. 1997 Numerical simulations of sink flow in the Hele-Shaw cell with small surface tension Eur. J. Appl, Math. 8, 533–550.
  * La Roche et al. (1991) La Roche, H., Fernández, J.F., Octavio, M., Loeser, A.G. & Lobb, C.J. 1991 Diffusion-limited-aggregation model for Poisson growth. Phys. Rev. A, 44, R6185(R).
  * Lakrout et al. (1999) Lakrout, H., Sergot, P. & Creton, C. 1999 Direct observation of cavitation and fibrillation in a probe tack experiment on model acrylic pressure-sensitive adhesives J. Adhesion, 69, 307–359.
  * Moffatt (1964) Moffatt, H.K. 1964 Viscous and resistive eddies near a sharp corner. J.Fluid Mech.  18, 1–18.
  * Moffatt (1977) Moffatt, H.K. 1977 Six lectures on general fluid dynamics and two on hydromagnetic dynamo theory. In Fluid Dynamics  (eds. R. Balian & J.-L. Peube), 149 – 233. Gordon and Breach.
  * Pihler-Puzović et al. (2012) Pihler-Puzović, D., Illien, P., Heil, M. & Juel, A. 2012 Suppression of complex fingerlike patterns at the interface between air and a viscous liquid by elastic membranes. Phys. Rev. Lett. , 108, 074502.
  * Poivet et al. (2004) Poivet, S., Nallet, F., Gay, C., Teissiere, J. & Fabre, P. 2004 Force response of a viscous liquid in a probe-tack geometry: fingering versus cavitation. Eur.Phys.J., E 15, 97–116.
  * Rabaut et al. (1988) Rabaut, M., Couder, Y. & Gerard, N. 1988 Dynamics and stability of anomalous Saffman-Taylor fingers. Phys. Rev. A 37, 935–947.
  * Saffman & Taylor (1958) Saffman, P.G. & Taylor, G.I. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. Lond.  A 245, 312–329.
  * Savina & Nepomnyashchy (2015) Savina, T.V. & Nepomnyashchy, A.A. 2015 On a Hele-Shaw flow with a time-dependent gap in the presence of surface tension. J. Phys. A: Math. Theor. 48, 125501 (13pp).
  * Savina (2021) Savina, T.V. 2021 On a two-fluid Hele-Shaw problem with an elliptical interface. J. Diff. Eqns. 270, 787–808.
  * Stefan (1874) Stefan, J. 1874 Veruche über die scheinbare Adhäsion. Sitzber. Akad. WisE. Wien, Math. naturw. Kl. A 69, 713.
  * Tanveer (2000) Tanveer, S. 2000 Surprises in viscous fingering. J.Fluid Mech. 409, 273-308.
  * Ward (2006) Ward, T. 2006 Radial spreading of a viscous drop between parallel-plane surfaces. Phys. Fluids  354, 816–824.
  * Ward (2011) Ward, T. 2011 Capillary-pressure driven adhesion of rigid-planar surfaces. J. Colloid Interface Sci.  18, 093101-8.
  * Webber & Huppert (2020) Webber, J.J. & Huppert, H.E. 2020 Time to approach similarity. Quart. J. Mech. Appl. Math. 73, 1-23.