# Compression-driven viscous fingering in a radial Hele-Shaw cell

Callum Cuttle  Liam C. Morrow  Christopher W. MacMinn
<EMAIL_ADDRESS>Department of Engineering Science,
University of Oxford, Oxford, OX1 3PJ, UK

###### Abstract

The displacement of a viscous liquid by a gas within a Hele-Shaw cell is a
classical problem. The gas-liquid interface is hydrodynamically unstable,
forming striking finger-like patterns that have been the subject of intense
focus. Generally, both the gas and liquid phases are taken to be
incompressible, with the capillary number being the key parameter that
determines the severity of the instability. Here, we consider a radially
outward displacement driven by the steady compression of a gas reservoir. The
associated gas injection rate is then unsteady due to the compressibility of
the gas. We identify a second nondimensional parameter, the compressibility
number, that plays a strong role in the development of the fingering pattern.
We use an axisymmetric model to study the impact of compressibility number on
the unsteady evolution of injection rate and gas pressure. We use linear
stability analysis to show that increasing the compressibility number delays
the onset of finger development relative to the corresponding incompressible
case. Finally, we present and compare a series of experiments and fully
nonlinear simulations over a broad range of capillary and compressibility
numbers. These results show that increasing the compressibility number
systematically decreases the severity of the fingering pattern at high
capillary number. Our results provide an unprecedented comparison of
experiments with simulations for viscous fingering, a comprehensive
understanding of the role of compressibility in unstable gas-liquid
displacement flows, and insight into a new mechanism for controlling the
development of fingering patterns.

##  I Introduction

When a fluid is displaced from a Hele-Shaw cell or porous medium by the
injection of a less viscous fluid, the fluid-fluid interface is
hydrodynamically unstable and tends to deform into complex, branched
structures [1, 2, 3, 4]. This classical viscous-fingering instability has been
extensively studied as an archetype of interfacial pattern formation [5, 6, 7]
and for its relevance to enhanced oil recovery [8], for which fingering poses
major obstacles. Modern applications include operation of fuel cells [9],
growth of bacterial colonies [10, 11], enhancing mixing in microfluidic
devices [12, 13], and dendrite solidification [14, 15]. Whether fingering is
advantageous or problematic, the key concern in all applications has been
understanding the mechanisms that influence the development of the fingering
pattern, which is driven by viscous forces in the defending phase and opposed
by capillary forces at the interface. The (modified) capillary number Ca [16],
which measures the relative scales of these forces, is therefore the key
control parameter. Recent studies have considered a variety of perturbations
to classical viscous fingering, such as imposing a time-dependent injection
rate [17, 18], replacing one of the rigid plates with an elastic membrane or
slab [19, 20], varying the gap between the plates as the gas is injected [21,
22], using a tapered (rather than uniform) flow cell [23, 24], employing
chemically reactive fluids with time-varying viscosities [25], and the
application of external electric fields [26]. Here, we examine a simpler but
thus far overlooked adaptation of the classical viscous fingering problem,
which arises when the injected fluid is a gas.

Gas-driven displacement of liquid is a common practical and experimental
scenario that is highly susceptible to viscous fingering. In addition to being
much less viscous than liquids, gases are also much more compressible than
liquids; as a result, they will undergo some amount of spring-like volumetric
compression under the typical viscous pressures of displacement flows. Such
compression-driven flows can exhibit unsteady flow rates [27, 28, 29], which
are known to exert a fundamental influence on pattern-forming processes [30].
Mathematical models of gas-liquid displacement in a Hele-Shaw cell typically
assume that the gas is incompressible, whereas experimentalists often take
pains to avoid or ignore gas compression, such as by withdrawing the liquid at
the outlet instead of injecting the air at the inlet [31, 32, 33, 34, 35, 36],
performing analysis based on the instantaneous interface velocity [37, 24], or
constraining themselves to low injection rates [38]. Alternatively, one may
inject a low-viscosity liquid rather than a gas, such as injecting water into
oil [39, 40]. However, the ability to entirely neglect the viscosity of the
gas phase makes gas-driven displacement flows especially tractable to
analytical and numerical analysis and many applications inherently involve
gases. As such, compressibility is typically considered to be a necessary
complication in studies of viscous fingering, but no study has yet offered
careful consideration of the impact of compressibility on displacement flows
and viscous fingering in Hele-Shaw cells.

Here, we use a combination of mathematical modelling, simulations, and
experiments to fully elucidate the role of gas compression during gas-driven
viscous fingering. We show that the unsteady flow rates observed in
experiments and simulations can be rationalised by a simple axisymmetric model
that couples the compression of an ideal gas to laminar viscous flow of
liquid. Within the context of this simplified model, the problem is controlled
by a single dimensionless compressibility number $\mathcal{C}$ that compares
the rate of viscous depressurisation to the rate of compressive
pressurisation. We find that $\mathcal{C}$ plays a similarly important role in
experiments and simulations involving fingering, dictating both the volumetric
growth rate and the evolution of pressure. We examine the dynamical-systems
framework that underpins the basic compression-driven displacement flow in the
simplified model, identifying the emergence of two distinct dynamical regimes
that are analogous to those recently described for compression-driven
displacement in a capillary tube [29]. A linear stability analysis of this
axisymmetric model suggests that compressibility can significantly delay the
onset of viscous fingering at high Ca, via the low initial injection rates of
the compression-driven flow. This prediction is confirmed in our experiments
and nonlinear simulations, where we observe a systematic delay in the onset of
fingering to larger radii as $\mathcal{C}$ is increased. We further show that
increasing $\mathcal{C}$ is as effective in delaying onset as decreasing Ca,
pointing to compressibility as a powerful new control parameter for gas-driven
viscous fingering. Our study also provides an unusually thorough comparison
between simulations and experiments of viscous fingering. Previous studies
have made quantitative comparisons between experimental results and
simulations over a relatively small portion of the parameter space [41, 17,
18, 42, 43]. The present study is the most extensive comparison of its type,
featuring quantitative examination of the fingering patterns, the time-
dependent injection rate, gas pressure, and the onset and development of the
instability. This work therefore provides a rare opportunity to validate a
numerical realisation of viscous fingering directly against a dedicated set of
experiments.

Figure 1: Schematic diagram of the fingering problem considered here. A
reservoir of gas (white) is compressed at a volumetric rate $Q$ in order to
displace liquid (blue) from a Hele-Shaw cell comprising the uniform gap $b$
between two plates.

##  II Mathematical Model

###  II.1 Full quasi-2D model

We consider a rigid, circular Hele–Shaw cell of radius $R$ and gap thickness
$b$, initially filled with an incompressible liquid of dynamic viscosity
$\mu$. Gas is injected through a hole in the centre of the cell by compressing
an attached gas reservoir of initial volume is $V_{\mathrm{res}}(0)$ at a
constant volumetric rate $Q$. We take the gas and liquid phases to be
immiscible. We ignore viscous pressure gradients within the gas phase, such
that the gas pressure $p_{g}(t)$ is spatially uniform.

Within the domain $\Omega$ of the liquid phase, we consider a quasi
two–dimensional (quasi-2D) model of Hele–Shaw flow that relates the gap-
averaged velocity of the liquid $\boldsymbol{v}$ to its pressure $p$ via

$\displaystyle\boldsymbol{v}$
$\displaystyle=-\frac{b^{2}}{12\mu}\mbox{\boldmath$\nabla$}p,\qquad$
$\displaystyle\boldsymbol{x}$
$\displaystyle\in\mathbb{R}^{2}\backslash\Omega(t).$ (1a)

Further, the incompressibility of the liquid phase requires that
$\mbox{\boldmath$\nabla$}\cdot\boldsymbol{v}=0$ and, hence, that

$\displaystyle\nabla^{2}p$ $\displaystyle=0,\qquad$
$\displaystyle\boldsymbol{x}$
$\displaystyle\in\mathbb{R}^{2}\backslash\Omega(t).$ (1b)

Following Peng  _et al._ [44], we impose the kinematic and dynamic conditions
at the moving interface $\partial\Omega(t)$:

$\displaystyle v_{n}$
$\displaystyle=-\frac{1}{1-f_{1}}\frac{b^{2}}{12\mu}\mbox{\boldmath$\nabla$}p\cdot\boldsymbol{n},\qquad\qquad$
$\displaystyle\boldsymbol{x}$ $\displaystyle\in\partial\Omega(t),$ (1c)
$\displaystyle p$
$\displaystyle=p_{g}(t)-\gamma\left(\frac{\pi}{4}\kappa+\frac{2f_{2}}{b}\right),\qquad\qquad$
$\displaystyle\boldsymbol{x}$ $\displaystyle\in\partial\Omega(t),$ (1d)

where $v_{n}$ is the normal velocity of the interface, $\boldsymbol{n}$ is the
unit normal to the interface, $\gamma$ is the interfacial tension, and
$\kappa$ is the signed in-plane curvature of the interface. The kinematic
boundary condition Eq. (1c) relates the velocity of the interface to the
velocity of the liquid. The dynamic boundary condition Eq. (1d) models the
capillary pressure via the Young–Laplace equation. Both of these boundary
conditions are influenced by the presence of thin residual films of liquid
that coat the walls in the gas region [45]. Following Peng  _et al._ [44], we
model these films via the empirical functions

$\displaystyle
f_{1}=\frac{|v_{n}\eta/\gamma|^{2/3}}{0.76+2.16|v_{n}\eta/\gamma|^{2/3}},\quad\mathrm{and}\quad
f_{2}=1+1.59|v_{n}\eta/\gamma|+\frac{|v_{n}\eta/\gamma|^{2/3}}{0.26+1.48|v_{n}\eta/\gamma|^{2/3}}.$
(1e)

We assume here and going forward that the liquid is perfectly wetting to the
cell walls. The total thickness of the films $f_{1}b$ depends on the normal
velocity of the interface at the instant of deposition, and we assume that the
films are static after deposition on the timescale of the flow, such that the
film thickness at each point in space is constant after the interface has
passed. We take the gas to be a fixed mass of ideal gas and the process to be
isothermal, meaning that all temperature changes due to compression and
expansion are rapidly equilibrated with the environment. Pressure and volume
are then related by Boyle’s law,

$\displaystyle p_{g}(t)V_{g}(t)$ $\displaystyle=p_{g}(0)V_{g}(0),$ (1f)

where $V_{g}(t)=V_{b}(t)+V_{\mathrm{res}}(t)$ is the total volume of gas in
the system, with $V_{b}$ the volume of gas in the cell (accounting for thin
films) and $V_{\mathrm{res}}$ the volume of the gas reservoir. Due to the
imposed steady compression rate $Q$, the volume of gas in the reservoir is
$V_{\mathrm{res}}(t)=V_{\mathrm{res}}(0)-Qt$, such that the pressure of the
gas is

$\displaystyle p_{g}(t)$
$\displaystyle=\frac{V_{g}(0)}{V_{\mathrm{res}}(0)-Qt+V_{b}(t)}\left[p_{\textrm{atm}}+\gamma\left(\frac{\pi}{4r_{0}}+\frac{2}{b}\right)\right],$
(1g)

where the term in square brackets is the initial gas pressure $p_{g}(0)$,
which balances the Laplace pressure jump across the interface via Eq. (1d) for
$p(0)=p_{\textrm{atm}}$. We have taken the system to be initially at rest,
with a small circular bubble of gas of radius $r_{0}$ centred on the origin
(see § II.3) such that the initial in-plane curvature of the interface is
$1/r_{0}$. At the outlet of the Hele–Shaw cell (around the rim), we take the
pressure of the liquid to be atmospheric,

$\displaystyle p$ $\displaystyle=p_{\textrm{atm}},\qquad|\boldsymbol{x}|\in
R.$ (1h)

The flow is driven by the eventual difference between the liquid pressure at
the liquid-gas interface and at the outlet.

We non-dimensionalise our problem via

$\displaystyle\boldsymbol{\hat{x}}=\frac{\boldsymbol{x}}{R},\quad\hat{t}=\frac{Q}{\pi
R^{2}b}t,\quad\hat{p}=\frac{\pi b^{3}}{6\mu
Q}p,\quad\hat{\boldsymbol{v}}=\frac{2\pi
Rb}{Q}\boldsymbol{v},\quad\hat{V}=\frac{V}{\pi R^{2}b},$ (2)

leading to

$\displaystyle\hat{\nabla}^{2}\hat{p}$ $\displaystyle=0,$
$\displaystyle\hat{\boldsymbol{x}}$
$\displaystyle\in\mathbb{R}^{2}\backslash\Omega(\hat{t}),$ (3a)
$\displaystyle\hat{v}_{n}$
$\displaystyle=-\frac{1}{1-f_{1}}\mbox{\boldmath$\nabla$}\hat{p}\cdot\boldsymbol{n},$
$\displaystyle\hat{\boldsymbol{x}}$ $\displaystyle\in\partial\Omega(\hat{t}),$
(3b) $\displaystyle\hat{p}$
$\displaystyle=\hat{p}_{g}(t)-\textrm{Ca}^{-1}\left(\frac{\pi}{4}\hat{\kappa}+2\alpha
f_{2}\right),\qquad\qquad$ $\displaystyle\hat{\boldsymbol{x}}$
$\displaystyle\in\partial\Omega(\hat{t}),$ (3c) $\displaystyle\hat{p}$
$\displaystyle=\mathcal{P},\qquad\qquad$ $\displaystyle|\hat{\boldsymbol{x}}|$
$\displaystyle\in 1,$ (3d) where $\displaystyle\hat{p}_{g}$
$\displaystyle=\frac{1}{1+\mathcal{V}^{-1}[\hat{V}_{b}(\hat{t})-\mathcal{S}^{2}-\hat{t}]}\left[\mathcal{P}+\textrm{Ca}^{-1}\left(\frac{\pi}{4\mathcal{S}}+2\alpha\right)\right]$
(3e) and $\mathcal{S}=r_{0}/R$ is the dimensionless initial radius of the
interface. The thin-film factors [Eq. (1e)] are then written as $\displaystyle
f_{1}=\frac{|\textrm{Ca}^{\prime}\hat{v}_{n}|^{2/3}}{0.76+2.16|\textrm{Ca}^{\prime}\hat{v}_{n}|^{2/3}},\quad\mathrm{and}\quad
f_{2}=1+1.59|\textrm{Ca}^{\prime}\hat{v}_{n}|+\frac{|\textrm{Ca}^{\prime}\hat{v}_{n}|^{2/3}}{0.26+1.48|\textrm{Ca}^{\prime}\hat{v}_{n}|^{2/3}}.$
(3f) The dimensionless control parameters are then
$\displaystyle\textrm{Ca}^{\prime}=\frac{\mu Q}{2\pi
Rb\gamma},\quad\mathcal{P}=\frac{\pi b^{3}p_{atm}}{6\mu
Q},\quad\mathcal{V}=\frac{V_{g}(0)}{\pi
R^{2}b},\qquad\alpha=\frac{R}{b},\qquad\mathrm{and}\qquad\mathcal{S}=\frac{r_{0}}{R},$
(3g)

where $\textrm{Ca}^{\prime}$ is the nominal capillary number and
$\textrm{Ca}=12\alpha^{2}\textrm{Ca}^{\prime}$ is the “modified” capillary
number, as used in previous work [e.g., 16]. From this point on, we refer to
Ca as the capillary number.

For comparison, we also consider the case where the gas is incompressible. In
this limit, Eq. (3e) is replaced by the constraint that $p_{g}$ must be chosen
so that the actual injection rate
$\hat{Q}_{b}\equiv\mathrm{d}\hat{V}_{b}/\mathrm{d}\hat{t}=1$ at all times and
the parameters $\mathcal{P}$ and $\mathcal{V}$ can be eliminated. This
constraint can be enforced by choosing $p_{g}$ such that

$\displaystyle\int_{0}^{2\pi}-\frac{\partial\hat{p}}{\partial\hat{r}}\bigg{\rvert}_{\hat{r}=1}\,\textrm{d}\theta=2\pi,$
(4)

which ensures that the total rate of outflow is unity.

###  II.2 Reduced model

To gain insight into the coupling between gas compression and liquid
displacement, we consider the axisymmetric limit in which the interface
remains circular with radius $r=s_{0}(t)$. Equations (3a)-(3e) then reduce to

$\displaystyle\frac{\text{d}\hat{s}_{0}}{\text{d}\hat{t}}$
$\displaystyle=-\frac{1}{2(1-f_{1})\hat{s}_{0}\ln\left(\hat{s}_{0}\right)}\left[\Delta\hat{p}_{g}-\textrm{Ca}^{-1}\left(\frac{\pi}{4\hat{s}_{0}}+2\alpha
f_{2}\right)\right],$ (5) $\displaystyle\Delta\hat{p}_{g}$
$\displaystyle=\frac{1}{1+\mathcal{V}^{-1}\left(\hat{s}_{0}(\hat{t})^{2}-\mathcal{S}^{2}-\hat{t}\right)}\left[\mathcal{P}+\textrm{Ca}^{-1}\left(\frac{\pi}{4\mathcal{S}}+2\alpha\right)\right]-\mathcal{P},$
(6)

where $\Delta\hat{p}_{g}=\hat{p}_{g}-\mathcal{P}$ is the dimensionless gauge
gas pressure. We can simplify further by taking $\mathcal{V}\gg 1$ and
$\mathcal{P}\gg\textrm{Ca}^{-1}[\pi/(4\mathcal{S})+2\alpha]$; the former limit
corresponds to an initial volume of gas that is much larger than the volume of
liquid in the cell (i.e., a large gas reservoir) while the latter corresponds
to negligible capillary pressure relative to atmospheric pressure, which also
implies $\mathcal{P}\gg\textrm{Ca}^{-1}[\pi/(4\hat{s}_{0})+2\alpha f_{2}]$.
Finally, ignoring thin film effects, such that $f_{1}=0$, Eqs. (5)-(6) reduce
to the ordinary differential equation

$\displaystyle\frac{\text{d}\hat{s}_{0}}{\text{d}\hat{t}}$
$\displaystyle=-\frac{\Delta\hat{p}_{g}}{2\hat{s}_{0}\ln\left(\hat{s}_{0}\right)},\quad\textrm{with}\quad\Delta\hat{p}_{g}=2\mathcal{C}^{-1}\left(\hat{t}-\hat{s}_{0}^{2}+\mathcal{S}^{2}\right).$
(7)

The evolution of the interface then depends on a single dimensionless
parameter

$\displaystyle\mathcal{C}=\frac{2\mathcal{V}}{\mathcal{P}}=\frac{12Q\mu
V_{g}(0)}{\pi^{2}R^{2}b^{4}p_{\textrm{atm}}},$ (8)

which we refer to as the compressibility number, and which is directly
analogous to the compressibility number identified in Ref. [29]. When the gas
is incompressible ($\mathcal{C}\to 0$), Eq. (7) degenerates to
$\hat{s}_{0}(\hat{t})=\left(\mathcal{S}^{2}+\hat{t}\right)^{1/2}$ and
$\Delta\hat{p}_{g}=-\ln(\hat{s}_{0})$, as expected.

Throughout the present work, we focus on injection via the steady compression
of a gas reservoir, as may be imposed by a syringe pump, for consistency with
our experiments. However, another common experimental implementation of gas
injection is from a high-pressure source regulated by a needle valve. In
Appendix A, we derive a mathematical model for this alternative approach and
show that the two methods are identical under the assumptions of the reduced
model.

###  II.3 Numerical scheme

We solve Eqs. (3a)-(3e) using the numerical scheme proposed by Morrow  _et
al._ [46], which we briefly summarize here. The scheme is based on the level-
set method [47], where we construct a level-set function $\phi$ such that
$\phi<0$ in the gas region and $\phi>0$ in the liquid region. We evolve $\phi$
via the level-set equation

$\displaystyle\frac{\partial\phi}{\partial\hat{t}}+F\left|\hat{\nabla}\phi\right|=0,$
(9)

where

$\displaystyle F=-\hat{\mbox{\boldmath$\nabla$}}\hat{p}\cdot\boldsymbol{n},$
(10)

and
$\boldsymbol{n}=\hat{\mbox{\boldmath$\nabla$}}\phi/|\hat{\mbox{\boldmath$\nabla$}}\phi|$
is the unit (outward) normal. This choice of speed function $F$ satisfies the
kinematic boundary condition on the interface [Eq. (3b)] and gives a
continuous expression for $F$ in the liquid region
$\boldsymbol{x}\in\mathbb{R}^{2}\backslash\Omega(\hat{t})$. We extend $F$ into
the gas region by solving the biharmonic equation [48]. To solve Eq. (9), we
use a second-order essentially non-oscillatory scheme for the spatial
derivatives, and integrate in time using second order total-variation-
diminishing Runge-Kutta with $\Delta t=\Delta x/[4\max|F|]$. We solve Eq. (3b)
for the liquid pressure via a finite difference stencil. Following Gibou  _et
al._ [49], we modify the stencil at nodes adjacent to the gas-liquid interface
to incorporate the dynamic boundary condition [Eq.(3c)], where the signed
curvature of the interface is
$\kappa=\mbox{\boldmath$\nabla$}\cdot\boldsymbol{n}$. Further, the volume of
the gas region is computed with

$\displaystyle\hat{V}_{g}=\int_{\mathbb{R}^{2}}(1-f)H(\phi)\,\textrm{d}\hat{V},$
(11)

where

$\displaystyle H=\begin{cases}0&\textrm{ if }\phi<-\delta\\\ 1&\textrm{ if
}\phi>\delta\\\
\frac{1}{2}\left(1+\frac{\phi}{\delta}+\frac{1}{\pi}\sin\frac{\pi\phi}{\delta}\right)&\text{
otherwise}\\\ \end{cases},$ (12)

$\delta=1.5\Delta x$, and $f$ is the proportion of the Hele-Shaw cell filled
with liquid at each node, as determined from $f_{1}$ [see Eq. (3f)]. When the
gas is incompressible, we discretise the integral in Eq. (4) via the
trapezoidal rule. We solve the resulting system of linear equations via LU
decomposition. All simulations are performed with the initial condition

$\displaystyle\hat{s}(\theta,0)=\mathcal{S}\left\\{1+\sum_{n=2}^{12}\varepsilon_{n}\cos\left[n(\theta-2\pi\theta_{n})\right]\right\\},$
(13)

where $\theta_{n}$ and $\varepsilon_{n}$ are selected at random from uniform
distributions on the intervals (0, 1) and ($5\times 10^{-4}$, $10^{-3}$),
respectively. Simulations are performed on the computational domain $0\leq
r\leq 1$ and $0\leq\theta<2\pi$ using $1000\times 1000$ equally spaced nodes.
Simulations are stopped when the maximum radius of the interface is 0.99; we
denote the time at which this occurs as $\hat{t}=\hat{t}_{f}$.

##  III Experimental methods

###  III.1 Set-up

Experiments are performed in a Hele-Shaw cell comprising two float-glass
plates of radius $R=105$ mm. We impose a gap $b=0.42\pm 0.01$ mm between the
plates using a plastic spacer. The spacer supports the outer 5 mm of the
plates and, in doing so, obstructs a small fraction of the outflow area. Note
that this partial obstruction is not included in the simulations, and is
believed to contribute to the systematically greater pressures recorded in
experiments [see Fig. 4(d)], although this discrepancy could instead result
from qualitative differences in fingering patterns (see § V.2 and § V.5).

The cell is initially filled with 1,000 cSt silicone oil (dynamic viscosity
$\mu=0.97$ Pa s and surface tension $\gamma=21$ mN m-1 at laboratory
temperature $22\pm 1^{\circ}$C; Sigma), which is dyed (Sudan III; Merck) and
filtered. We impose a fixed hydrostatic pressure at the outer rim of the cell
(the outlet) by surrounding the cell with a shallow oil reservoir, filled to a
depth approximately 1 mm above the top of the gap, which is maintained using
an overflow.

The cell is connected to an air reservoir via a 2 mm diameter injection port.
The reservoir comprises two 50 mL airtight glass syringes (1050TTL; Hamilton),
along with stiff connective tubing of internal volume $12\pm 1$ mL (Legris).
Short sections of connective tubing (Tygon) are used sparingly to minimise
pressure-induced changes in the volume of the air reservoir itself. The total
initial volume of the air reservoir (tubing and syringes) is set to
$V_{\mathrm{res}}(0)\in\\{25,50,100,200\\}$ mL, with a relative error of 2-4%.
To achieve $V_{\mathrm{res}}(0)$=200 mL, we connect an additional sealed
acrylic box of internal volume $120\pm 3$ mL. The total internal volume of the
box and tubing was measured via changes in air pressure during controlled
compression tests with the system closed. These tests also suggest that air
leakage was negligible over the timescales and pressures of our experiments.

Prior to each experiment, we introduced a circular precursor bubble of initial
radius $s_{0}(0)=2.7\pm 0.1$ mm by injecting air very slowly, such that no
significant pre-compression of the air was introduced. The initial air
pressure is taken to be atmospheric, neglecting the small hydrostatic pressure
imposed by the oil reservoir and the Laplace pressure jump at the interface.
The initial volume of the bubble, around 0.01 mL, is negligible compared with
$V_{\mathrm{res}}(0)$.

To conduct an experiment, the air was then compressed using a syringe pump
(AL-4000; WPI) at a fixed volumetric rate $Q\in\\{1.25,2.50,5.00,10.0\\}$ mL
min-1. The gauge pressure $\Delta p_{g}$ of the air relative to atmospheric
pressure $p_{\textrm{atm}}$ was recorded using a pressure sensor (0-15 PSI;
Honeywell) via a USB DAQ (U6; LabJack) at a frequency of approximately 20 Hz.
The pressurised air drove oil out of the cell and into the surrounding oil
reservoir.

We imaged the motion of the interface using a CMOS camera (acA4096-30um;
Basler) mounted vertically below the cell and recording at a fixed frame rate
of 1.5–12 frames per second (fps), depending on $Q$, and at a spatial
resolution 91 $\mu$m per pixel. Over each experiment, we typically recorded
400 frames. The cell was backlit using a custom array of LEDs (Wholesale LED
Lights), diffused through opalescent acrylic (Sheet Plastics) and a blue
filter (Stage Depot) to enhance contrast with the dyed oil. Ambient light from
the laboratory was blocked by cloaking the set-up in opaque fabric (BK5;
Thorlabs).

Each experiment was performed twice to assess reproducibility. We found that
bulk displacement measurements, such as the evolution of the injection rate
and the air pressure, were highly reproducible between experiments despite
significant variation in the fingering patterns and their associated metrics
after onset (see § V.5 and Appendix B).

###  III.2 Data processing

The volume of the air in the cell $V_{b}$ was calculated from the air
pressure. We did so by modelling the air as a fixed mass of isothermal ideal
gas, such that $[p_{\textrm{atm}}+\Delta
p_{g}(t)]V_{g}(t)=p_{\textrm{atm}}V_{\mathrm{res}}(0)$. Taking
$V_{g}(0)=V_{\textrm{res}}(0)$ and $p_{g}(0)=p_{\textrm{atm}}$ by neglecting
the initial volume $V_{b}(0)$ of the bubble and the initial Laplace and
hydrostatic pressure, respectively, introduces negligible error in this
calculation. The total volume of air, which changes due to both compression of
the reservoir and invasion into the flow cell is
$V_{g}(t)=V_{\mathrm{res}}(0)-Qt+V_{b}(t)-V_{b}(0)$, so that

$\displaystyle V_{b}(t)=V_{b}(0)+Qt-\left[\frac{V_{\mathrm{res}}(0)\Delta
p_{g}(t)}{p_{\textrm{atm}}+\Delta p_{g}(t)}\right].$ (14)

Accounting for $V_{b}(0)$ at this step ensures that the initial value of
$V_{b}(t)$ is accurate. The actual time-dependent injection rate of air into
the cell $Q_{b}=\mathrm{d}V_{b}/\mathrm{d}t$ (which is equivalent to the flow
rate of liquid out of the cell) was calculated at each recorded $V_{b}(t)$
data point by taking a second-order polynomial least-squares fit to the data
on either side of that point. The first derivative of the fitted function was
then taken as the local value of $Q_{b}(t)$; the size of the fitting window
was automatically adjusted to minimise a $\chi^{2}$ fitting parameter that
avoided over- or underfitting to the data.

Recorded experimental images were processed in MATLAB R2020b to recover the
two-dimensional area of the bubble and perimeter of the interface as functions
of time. A reference image of the cell prior to air injection was subtracted
from every subsequent image to facilitate isolation of the interface. Pixels
along the interface were identified using the edge function with a canny
filter. The resulting binary image was dilated to ensure connectivity of the
interface, which was taken to be the largest connected region. A binary image
of the air region was generated by filling the connected interface contour and
then eroding to compensate for the earlier dilation step. The area and
perimeter of the air region were then calculated using regionprops and
converted to dimensional units according to the spatial resolution of the
camera. We process images up to the first frame when the maximum radial
coordinate of the interface $\max(s)$ is greater than $0.9R$, beyond which
interaction between the interface and the spacer became visually noticeable.
We refer to the moment when this occurs as the near-breakout time $t_{0.9}$.

##  IV Dimensional and non-dimensional parameter ranges

Going forward, we discuss results primarily in terms of Ca and $\mathcal{C}$.
In both experiments and simulations, however, we choose to vary the
dimensional nominal flow rate $Q$ and initial gas volume $V_{g}(0)$ for
practical reasons. Moreover, $\mathcal{C}$ does not appear directly in the
full model [Eqs. (3a)-(3e)], but is rather expressed as
$\mathcal{C}=2\mathcal{V}/\mathcal{P}\propto QV_{g}(0)$. To vary $\mathcal{C}$
while keeping $\textrm{Ca}\propto Q$ fixed, we fix $Q$ and vary the initial
gas volume within the range $V_{g}(0)\in${3.125, 6.25, 12.5, 25, 50, 100, 200,
400, 800, 1600} mL, equivalent to $\mathcal{V}\in${0.215, 0.430, 0.859, 1.72,
3.44, 6.87, 13.7, 27.5, 55, 110} (experiments were only performed at bold
values). To vary Ca while keeping $\mathcal{C}$ fixed, we vary $Q$ while
fixing the product $QV_{g}(0)$ (or, equivalently, $\mathcal{V}/\mathcal{P}$).
We restrict our simulations to the experimental values of $Q$, corresponding
to $\textrm{Ca}\in\\{2.61\times 10^{3}$, $5.21\times 10^{3}$, $1.04\times
10^{4}$, $2.08\times 10^{4}$}. Values of $\mathcal{P}$ and $\mathcal{V}$ are
given in the caption of each figure, where appropriate. The remaining
parameters are fixed at $\alpha=250$ and $\mathcal{S}=0.025$, with dimensional
values listed in § III.1.

##  V Results

Section V is organised as follows. We begin in § V.1 by describing the key
observations of compression-driven viscous fingering. In § V.2, we then
consider the ‘bulk’ displacement dynamics, and specifically the unsteady flow
rate and injection pressure, that arise from compression-driven displacement
in a radial Hele-Shaw cell. We show that the reduced model derived in § II.2
qualitatively captures the variation in bulk displacement dynamics with
increasing $\mathcal{C}$. Comparison between experiments, simulations and the
reduced model demonstrates that the bulk displacement dynamics of the full
system are also controlled primarily by $\mathcal{C}$ and are remarkably
insensitive to both variations in Ca and the presence of viscous fingering. In
§ V.3, we examine the underlying dynamical-systems structure that dictates the
unsteady flow rate in the reduced model. The dynamical regimes described in
this framework qualitatively predict the bursts of high flux observed at
$\mathcal{C}\gtrsim 1$ both in experiments and simulations. In § V.4, we
perform a linear stability analysis. At sufficiently large Ca, our analysis
suggests that compressibility significantly delays the onset of viscous
fingering relative to an incompressible system; we also identify a weaker
effect at low Ca, where compressibility can instead promote onset to
comparatively smaller radii. In § V.5, we examine the growth of viscous
fingers in both experiments and simulations. We observe a strong delay in the
onset with increasing $\mathcal{C}$, consistent with the predictions of our
linear stability analysis. Finally, we discuss the qualitative impacts of
compressibility on the nonlinear evolution of the fingering pattern in both
experiments and simulations.

###  V.1 Features of compression-driven viscous fingering

Figure 2: Experimental fingering patterns for Ca increasing top to bottom and
$\mathcal{C}$ increasing left to right. The corresponding dimensional
parameters are (top to bottom) $Q=[1.25,2.5,5,10]$ mL/min and (along each
north-west/south-east diagonal) $V_{\mathrm{res}}(0)=[25,50,100,200]$ mL, with
$R=10.5$ cm, $\mu=0.97$ Pa s, $\gamma=21$ mN m-1, and $b=0.42$ mm.

We begin by briefly examining the qualitative impact of compressibility on
viscous fingering. In Figs. 2 and 3, we show the evolution of the fingering
pattern across all of our experiments and simulations, respectively. For the
experiments [Fig. 2], each image includes the observed interface at time
intervals of $t_{0.9}/10$, where $t_{0.9}$ is the near-breakout time at which
the experiment is concluded (§ III.2). Rows correspond to fixed Ca, while
columns correspond to fixed compressibility number. In terms of dimensional
quantities, we can interpret this arrangement as each row having a fixed
nominal injection rate $Q$ with the initial gas reservoir size increasing left
to right. The behaviour for increasing Ca (top to bottom) is as expected from
classical work on viscous fingering [39]. The fingering pattern becomes more
severe as Ca increases in the sense that we observe more and narrower fingers,
as well as increasing instances of tip-splitting and side branching.
Additionally, the onset of fingering (i.e., the point at which the interface
deviates noticeably from a circle) appears to occur at smaller radii for
larger Ca. Both of these features are well established in previous work.

It is striking, however, that we see a similar decrease in onset radius with
decreasing $\mathcal{C}$ at fixed Ca (i.e., decreasing the air reservoir
volume while fixing the nominal injection rate), which corresponds to moving
from right to left along a given row. In other words, increasing $\mathcal{C}$
at fixed Ca, in this case by using a larger air reservoir, appears to
systematically delay the onset of viscous fingering. Varying the initial
reservoir volume $V_{\mathrm{res}}(0)$ does not change Ca, so its impact on
viscous fingering is not considered in classical studies. Rather, changing
$V_{\textrm{res}}(0)$ changes how ‘compressible’ the system is [30], as
measured by the value of $\mathcal{C}$. This change in compressibility
modifies the actual time-dependent injection rate (as distinct from the
nominal injection rate $Q$) due to coupling between viscous forces in the
displaced liquid and compressive forces in the air [29]. Our focus below is to
formally rationalise and quantify the impact of compressibility in this
system.

In addition to our experimental results, we conducted an extensive set of
simulations over a wider range of parameters. These simulation results are
shown in Figure 3, again with rows and columns ordered by Ca and
$\mathcal{C}$, respectively. The first column shows the results for an
incompressible system. Images with a shaded square background indicate
simulations conducted at the same parameters as the experiments in Fig. 2. Our
simulations and experiments are in strong qualitative agreement, showing
consistent variations in patterns as Ca and $\mathcal{C}$ are varied. We
analyse and compare these results quantitatively in § V.5.

Figure 3: Example numerical solutions for Ca increasing top to bottom
($\mathcal{P}=[194,96.9,48.5,24.2]$). Column one shows solutions of the
incompressible model, while columns two to eight show numerical solutions of
the compressible model Eqs. (3a)-(3e) with $\mathcal{C}$ increasing left to
right. Profiles are plotted in equal time intervals of $t_{f}/10$. Simulations
with shaded background correspond to the same parameter values as the
experiments shown in Fig. 2.

###  V.2 Bulk displacement dynamics

Figure 4: Compressibility number $\mathcal{C}$ has a strong impact on the
evolution of injection rate $\hat{Q}_{b}$ and gas gauge pressure
$\Delta\hat{p}_{g}$. Here, we plot (a) $\hat{Q}_{b}$ and (c)
$\Delta\hat{p}_{g}$ from from the numerical solution to the reduced model [Eq.
(7)] with $\mathcal{S}=0.025$ and (blue to red) $\mathcal{C}=0.018$, $0.036$,
$0.071$, $0.14$, $0.28$, $0.57$, $1.13$, $2.27$, and $4.54$. The
incompressible solution is shown for comparison (dashed black). Panels (b) and
(d) show the same quantities from numerical simulations (solid curves) and
experiments (symbols) with viscous fingering at $\textrm{Ca}=2.08\times
10^{4}$ ($Q=10$ mL/min and $V_{\mathrm{res}}(0)\in\\{25,50,100,200\\}$ mL; row
4 of Figs. 2 and 3). Blue and red shades indicate $\mathcal{C}<1$ and
$\mathcal{C}\geq 1$, respectively.

The dynamics of the reduced model [Eq. (7)] are illustrated in Fig. 4, which
shows the evolution of the nondimensional injection rate $\hat{Q}_{b}$ [Fig.
4(a)] and the gauge gas pressure $\Delta\hat{p}_{g}$ [Fig. 4(c)] for different
compressibility numbers $\mathcal{C}$. For reference, the incompressible
solution for a circular interface is also shown (dashed black), for which
$\hat{Q}_{b}=1$ and $\Delta\hat{p}_{g}$ monotonically decreases as liquid
drains and viscous resistance decreases. In the reduced model, by contrast,
the injection rate varies strongly as the interface advances and the gas
pressure evolves non-monotonically. This behaviour is due to the basic
coupling between viscous displacement and compressive pressurisation
(fingering is absent from this model). Initially, the gas compresses and
pressurises, such that $\hat{Q}_{b}$ and $\Delta\hat{p}_{g}$ increase
gradually from initial values of zero. As the interface advances and drives
liquid out, the resistance to flow decreases. The injection rate eventually
exceeds the nominal injection rate ($\hat{Q}_{b}>1$), at which point the gas
begins to expand and $\Delta\hat{p}_{g}$ decreases. For $\mathcal{C}\ll 1$,
the compressible dynamics differ only weakly from the incompressible dynamics,
with $\hat{Q}_{b}$ rapidly reaching and then exceeding the nominal flux,
before relaxing back toward the incompressible solution $\hat{Q}_{b}=1$. As
$\mathcal{C}$ approaches unity, the maximum in $\hat{Q}_{b}$ increases and
occurs at later times. For $\mathcal{C}>1$, the maximum in $\hat{Q}_{b}$
vanishes, and the injection rate instead increases monotonically, diverging as
the interface escapes the cell at the moment of breakout (when
$\hat{s}_{0}=1$). The qualitative change in dynamics around $\mathcal{C}=1$ is
consistent with the distinct dynamical regimes observed for compression-driven
displacement in a capillary tube [29], as discussed further in § V.3.

We show in Figure 4(b, d) that the full numerical simulations (solid lines)
and the experiments (symbols) with viscous fingering both exhibit
qualitatively similar displacement dynamics in $\hat{Q}_{b}$ and
$\Delta\hat{p}_{g}$. The presence of fingers in the experiments and
simulations leads to an earlier breakout than in the reduced model. In
addition, the transition from non-monotonic to monotonically increasing
$\hat{Q}_{b}$ occurs at a slightly lower value of $\mathcal{C}$ in the
simulations and experiments than in the reduced model. Otherwise, the
evolution of $Q_{b}$ and $\Delta\hat{p}_{g}$ and the variation with
$\mathcal{C}$ are strikingly similar. This agreement suggests that viscous
fingering has only a weak influence on the underlying displacement dynamics,
primarily leading to earlier breakout compared with the reduced model [Eq.
(7)]. Furthermore, the quantitative agreement between experiments and
corresponding simulations is in contrast with the visually distinct fingering
patterns generated in each case; comparing the bottom rows of Figs. 2 and 3,
we see that experiments produce fingering patterns with greater instances of
tip-splitting and side-branching, resulting in more severely distorted
interfaces (see § V.5).

Figure 5: Capillary number Ca has a modest impact on the evolution of
injection rate $\hat{Q}_{b}$ and gas gauge pressure $\Delta\hat{p}_{g}$. Here,
we plot (a) $\hat{Q}_{b}$ and (b) $\Delta\hat{p}_{g}$ for $\mathcal{C}=0.14$
from numerical simulations of viscous fingering at
$\textrm{Ca}\in\\{2.61\times 10^{3}$, $5.21\times 10^{3}$, $1.04\times
10^{4}$, $2.08\times 10^{4}$} ($\mathcal{P}\in\\{194,96.9,48.5,24.2\\}$,
$\mathcal{V}\in\\{1.72,3.44,6.87,13.7\\}$; column 5 in Fig. 3) and from the
reduced model (dashed black).

As shown in Fig. 5, the evolution of $\hat{Q}_{b}$ [Fig. 5(a)] and
$\Delta\hat{p}_{g}$ [Fig. 5(b)] at fixed $\mathcal{C}$ is also only weakly
modified by varying Ca over an order of magnitude, despite the strong
variation in fingering patterns (see fourth column in Fig. 2 and fifth column
of Fig. 3). The dominant effect of fingering on compression-driven
displacement is to induce an increasingly early breakout as Ca increases.
Larger values of Ca also lead to systematically lower values of
$\Delta\hat{p}_{g}$, consistent with the fact that more severe fingering
patterns bypass an increasingly large fraction of the liquid. Hence, the bulk
displacement dynamics, in terms of injection rates and driving pressures, are
remarkably insensitive to viscous fingering, and are governed primarily by
$\mathcal{C}$. As noted in § III, the systematically greater pressure recorded
in experiments [Fig. 4(d)] may derive in part from the geometry of the spacer
used to impose the gap, which was not accounted for in the simulations.
However, differences in the fingering patterns [Figs. 2 and 3] may also
influence the pressure evolution, as illustrated by the results of Fig. 5(b)
(see § V.5).

###  V.3 Dynamical systems framework

The reduced model exhibits two dynamical regimes, characterised by the change
from non-monotonic to monotonically increasing $\hat{Q}_{b}$ around
$\mathcal{C}\approx 1$ [Fig. 4$(a)$]. We next rationalise these regimes in
terms of a general dynamical systems framework. A similar treatment was
originally applied to compression-driven displacement in a capillary tube by
Cuttle and MacMinn [29]; we take the same approach here to describe the
dynamics of axisymmetric displacement in a Hele-Shaw cell, as embodied by the
reduced model [Eq. (7)].

We start by considering the injection rate, which for a circular interface is
$\hat{Q}_{b}=2\hat{s}_{0}(\mathrm{d}\hat{s}_{0}/\mathrm{d}\hat{t})$. From the
reduced model, $\hat{Q}_{b}=\Delta\hat{p}_{g}/\hat{\omega}$, where we have
introduced the resistance $\hat{\omega}=\ln\left(1/\hat{s}_{0}\right)$ by
analogy with Ohm’s law. We can then write the reduced model as

$\displaystyle\mathcal{R}(\hat{Q}_{b})=\mathcal{R}(\Delta\hat{p}_{g})-\mathcal{R}(\hat{\omega})=\frac{2}{\hat{\omega}\mathcal{C}}\left(\frac{1}{\hat{Q}_{b}}-1+\frac{\mathcal{C}\hat{Q}_{b}}{4\hat{s}_{0}^{2}}\right),$
(15)

where $\mathcal{R}(\hat{x})=\dot{x}/\hat{x}$ is the relative rate of change of
the variable $\hat{x}$, with $\dot{x}=\mathrm{d}\hat{x}/\mathrm{d}\hat{t}$.
The reduced model admits two trivial solutions, for which the driving
compressive force and the opposing viscous resistance decrease at the same
relative rate. These are

$\displaystyle\bar{Q}_{\pm}=2\hat{s}_{0}\left(\frac{\hat{s}_{0}\pm\sqrt{\hat{s}_{0}^{2}-\mathcal{C}}}{\mathcal{C}}\right),$
(16)

which satisfy $\mathcal{R}(\hat{Q}_{b})=0$.

Figure 6: Phase-space representation of the reduced model [Eq. (7)], in terms
of injection rate $\hat{Q}_{b}$ and interface radius $\hat{s}_{0}$, for
$\mathcal{C}=0.25$. The local time derivative of flow rate $\dot{Q}_{b}$ is
indicated by the colormap (see colorbar). Stable and unstable trivial
solutions $\bar{Q}_{\pm}$ are shown as solid and dashed black lines,
respectively. The $\hat{s}_{0}$ and $\hat{Q}_{b}$ coordinates of the saddle-
node bifurcation point are indicated by dot-dashed and dotted lines,
respectively. A solution of the reduced model for $\mathcal{C}=0.25$ and
$\mathcal{S}=0.01$ is plotted as a thick blue curve. Regions of the phase
space for which $\dot{Q}_{b}$ is positive or negative are shaded orange or
purple, respectively (see colorbar). When
$\mathcal{R}(\Delta\hat{p}_{g})>\mathcal{R}(\hat{\omega})$, the interface
accelerates, $\dot{Q}_{b}>0$. This occurs in three distinct subregions (see
inset): when $\hat{Q}_{b}<1$ because the gas is compressing while resistance
decreases ($\Delta\dot{p}_{g}>0$ and $\dot{\omega}<0$); when
$1<\hat{Q}_{b}<\bar{Q}_{-}$ because the gas is decompressing slowly
($\Delta\dot{p}_{g}<0$ is small); and when $\hat{Q}_{b}>\bar{Q}_{+}$ because
the gas is over-compressed ($\Delta\hat{p}_{g}$ is large). When
$\mathcal{R}(\Delta\hat{p}_{g})<\mathcal{R}(\hat{\omega})$, the interface
decelerates, $\dot{Q}_{b}<0$. This occurs only in the region
$\bar{Q}_{-}<\hat{Q}_{b}<\bar{Q}_{+}$, where the gas decompresses quickly
($\Delta\dot{p}_{g}<0$ and $|\Delta\dot{p}_{g}|$ is large).

The influence of these trivial solutions on compressible displacement dynamics
is illustrated by the phase-space plot in Fig. 6, where
$\bar{Q}_{\pm}(\hat{s}_{0})$ are plotted as dashed and solid black curves,
respectively, for $\mathcal{C}=0.25$. The trivial solutions separate regions
where $\dot{Q}_{b}$ is positive (orange) or negative (purple). The local
change of sign in $\dot{Q}_{b}$ determines the stability of each trivial
solution. Specifically, $\bar{Q}_{+}$ is a repeller and $\bar{Q}_{-}$ is an
attractor, in the sense that small perturbations grow or decay, respectively.
A solution of the reduced model is plotted as a faint blue curve, displaying
the characteristic nonmonotonic variation in $\hat{Q}_{b}$ with $\hat{s}_{0}$
and showing that, rather than tending back toward the incompressible solution
($\hat{Q}_{b}=1$), the system is drawn onto the attractive solution
$\bar{Q}_{-}$.

Due to their $\hat{s}_{0}$-dependence, the trivial solutions only exist for
$\hat{s}_{0}\geq\sqrt{\mathcal{C}}$. At $\hat{s}_{0}=\sqrt{\mathcal{C}}$, the
two branches $\bar{Q}_{\pm}$ meet and annihilate at a saddle-node bifurcation.
For $\hat{s}_{0}<\sqrt{\mathcal{C}}$ there are no trivial solutions:
$\dot{Q}_{b}>0$ for all $\hat{Q}_{b}$, so that the interface accelerates
monotonically. Because $\hat{s}_{0}$ must be less than unity, $\mathcal{C}>1$
implies that there is no attractive solution and the interface accelerates
monotonically toward breakout. Even for $\mathcal{C}<1$, however, the dynamics
may fail to converge onto $\bar{Q}_{-}$, depending on the initial radius
$\mathcal{S}$. Computing the basin of attraction in terms of $\mathcal{S}$ for
a given $\mathcal{C}$ is only possible numerically, and is beyond the scope of
this study. A brief exploration suggests that, for $\mathcal{S}=0.01$, the
system fails to converge onto $\bar{Q}_{-}$ when $\mathcal{C}\geq 0.92$.
Therefore, $\mathcal{C}=1$ is an upper bound on the critical value at which
the stability of the flow changes, unlike in a capillary tube where the
critical value $\mathcal{C}=1$ is a precise indicator of the flow regime [29].

The stable trajectories that do converge onto $\bar{Q}_{-}$ terminate with a
breakout flux of $\hat{Q}_{b}(\hat{s}_{0}=1)\equiv
2\left(1-\sqrt{1-\mathcal{C}}\right)/\mathcal{C}$. This is because the
resistance $\hat{\omega}$ vanishes, driving the system exactly onto
$\bar{Q}_{-}$ at the moment of breakout. Similarly, for unstable trajectories
that fail to converge onto $\bar{Q}_{-}$, the vanishing resistance drives
divergent $\mathcal{R}(\hat{Q}_{b})$, and hence divergent $\hat{Q}_{b}$, in
the absence of a local attractive solution.

Figure 7: (a) Breakout time $\hat{t}_{f}$ at which the interface first reaches
the rim of the flow cell and (b) corresponding breakout gauge gas pressure
$\Delta\hat{p}_{g}(\hat{t}_{f})$, both computed from numerical solution of Eq.
(7) for $\mathcal{S}=0.025$ (black lines) and from the full numerical
simulations (colored curves and symbols) for $\textrm{Ca}\in\\{2.61\times
10^{3}$, $5.21\times 10^{3}$, $1.04\times 10^{4}$, $2.08\times 10^{4}$}
(increasing in the direction of the arrows). Each point is the average value
from 10 simulations, with error bars equal to one standard deviation above and
below the mean. Error bars are smaller than symbols in panel (a).

In Figures 7(a, b), we plot the breakout time and the breakout pressure,
respectively. We compute these values numerically from the reduced model
(solid lines) for $\mathcal{S}=0.025$. We find that, for $\mathcal{C}\lesssim
1$, the breakout time $\hat{t}_{f}$ is almost exactly 1 (within numerical
resolution) while the breakout pressure $\Delta\hat{p}_{g}(\hat{t}_{f})$ is
almost exactly zero. This corresponds to a scenario where breakout occurs at
exactly the same time as for an incompressible flow driven at the nominal flow
rate. As a consequence, the volume of air displaced by the piston at the
moment of breakout is exactly equal to the volume of liquid displaced by the
air, such that the air returns to its initial volume and pressure. Hence, the
driving compressive force vanishes at the same rate as the viscous resistive
force, consistent with terminating on the trivial solution $\bar{Q}_{-}$. For
$\mathcal{C}\gtrsim 1$, in contrast, breakout is delayed ($\hat{t}_{f}>1$)
which means that the air is still compressed at breakout and
$\Delta\hat{p}_{g}(\hat{t}_{f})>0$. Hence, the driving pressure remains finite
as the opposing resistance vanishes and $\hat{Q}_{b}$ diverges.

For comparison, in Fig. 7 we also show the breakout time and the breakout
pressure from simulations with fingering over a range of $\mathcal{C}$ at
different values of Ca (colored curves and symbols). We observe qualitatively
similar behavior, with $\hat{t}_{f}$ and $\Delta\hat{p}_{g}(\hat{t}_{f})$
varying slowly for $\mathcal{C}\lesssim 0.1$ before increasing sharply around
$\mathcal{C}\approx 0.1$. We also observe nonmonotonic variation in
$\Delta\hat{p}_{g}(\hat{t}_{f})$ for the highest Ca studied, and we would
expect to see the same for the lowest Ca were the range of $\mathcal{C}$
extended. Breakout times are systematically and substantially lower in the
simulations than in the reduced model because a significant fraction of the
liquid is bypassed by viscous fingers, allowing the interface to reach the
edge of the cell earlier. Similarly, the breakout pressures are systematically
greater in the fingering simulations, which may reflect the significant volume
of liquid left in the cell at the moment of breakout; because the interface is
still advancing, a significant pressure is still required to drive flow in the
remaining liquid.

###  V.4 Linear stability analysis

Figure 8: (a) Growth rate of the most unstable mode of perturbation [Eq. (19)]
and (c) most unstable mode of perturbation [Eq. (20)] for $\mathcal{C}=0.14$
and $\textrm{Ca}\in\\{2.61\times 10^{3}$, $5.21\times 10^{3}$, $1.04\times
10^{4}$, $2.08\times 10^{4}\\}$ (increasing in the direction of the arrows).
Panels (b) and (d) show the same quantities for $\textrm{Ca}=1.04\times
10^{4}$ and (blue to red) $\mathcal{C}\in\\{0.018$, $0.036$, $0.071$, $0.14$,
$0.28$, $0.57$, $1.13$, $2.27\\}$ (increasing in the direction of the arrows).
Blue and red shades indicate $\mathcal{C}<1$ and $\mathcal{C}\geq 1$,
respectively, and the classical incompressible values are shown for comparison
(dashed black). Note that while the incompressible values are defined for all
$\hat{s}_{0}\geq 0$, compressible values are only defined for
$\hat{s}_{0}\geq\mathcal{S}=0.025$.

We next examine the impact of compressibility on the onset of viscous
fingering by performing a linear stability analysis. To do so, we consider a
slightly perturbed circular solution of the form

$\displaystyle\hat{s}(\theta,\hat{t})$
$\displaystyle=\hat{s}_{0}(\hat{t})+\varepsilon\hat{\gamma}_{n}(\hat{t})\cos
n\theta+\mathcal{O}(\varepsilon^{2}),$ (17)
$\displaystyle\hat{p}(\hat{r},\theta,\hat{t})$
$\displaystyle=\hat{p}_{0}(\hat{r},\hat{t})+\varepsilon\hat{A}_{n}(\hat{r},\hat{t})\cos
n\theta+\mathcal{O}(\varepsilon^{2}),$ (18)

where $\varepsilon\ll 1$, $n\geq 2$, and $\hat{s}_{0}$ and $\hat{p}_{0}$ are
the unperturbed circular solution (i.e., the base state). Here,
$\hat{\gamma}_{n}$ and $\hat{A}_{n}$ denote, respectively, the amplitudes of
the $n$th mode of perturbation to the radius and pressure. Following Paterson
[5], we use the $\mathcal{O}(\varepsilon)$ problem to derive an evolution
equation for the relative growth rate of the $n$th mode of perturbation,

$\displaystyle\mathcal{R}(\hat{\gamma}_{n})=\frac{n-1}{2\hat{s}_{0}^{2}}\left[\hat{Q}_{b}-\frac{n(n+1)}{\textrm{Ca}\hat{s}_{0}}\right].$
(19)

The most unstable mode of perturbation is then

$\displaystyle n_{\max}=\sqrt{\frac{1+\textrm{Ca}\hat{s}_{0}\hat{Q}_{b}}{3}},$
(20)

which comes about by solving $\partial\mathcal{R}(\hat{\gamma}_{n})/\partial
n=0$. Thus, the value of $n_{\max}$ depends on the evolution of
$\hat{Q}_{b}(\hat{s}_{0})$ for an unperturbed circular interface. This base
state is precisely the solution to the reduced model, where
$\hat{Q}_{b}=2\hat{s}_{0}\dot{s}_{0}$. By combining the reduced model and
linear stability analysis, we may therefore understand how compression-driven
displacement modifies the onset of viscous fingering.

Figures 8(a, c) show $\mathcal{R}(\hat{\gamma}_{\max})$ and $n_{\max}$,
respectively, as functions of $\hat{s}_{0}$ for fixed $\mathcal{C}$ at
different values of Ca. These curves are calculated by solving Eq. (7)
numerically and substituting $\hat{Q}_{b}(\hat{s}_{0})$ into Eqs. (19) and
(20). Increasing Ca enhances the growth rate of the instability
[$\mathcal{R}(\hat{\gamma}_{n})$ increases] as well as the most unstable mode
of perturbation $n_{\max}$ for all $\hat{s}_{0}$. These observations are
consistent with our experimental and numerical results (Figs. 2 and 3) where,
over the values of $\mathcal{C}$ considered, increasing Ca results in more
prominent branching and tip splitting behaviour.

Fixing Ca and vary $\mathcal{C}$, as in Figs. 8(b, d), we find that increasing
$\mathcal{C}$ suppresses the instability at earlier times and promotes it at
later times. This observation is consistent with the fact that, as discussed
in section V.2, compression-driven displacement is characterised by a lower
injection rate at earlier times and a greater injection rate at later times,
compared with the nominal injection rate. As a result, for a given
$\mathcal{C}$, both $\mathcal{R}(\hat{\gamma}_{n})$ and $n_{\max}$ cross the
incompressible solution (dashed line; $\hat{Q}_{b}=1$) at some value of
$\hat{s}_{0}$. Before this crossover, the instability is suppressed relative
to an incompressible flow at the same Ca; after this crossover, the
instability is promoted.

Physically, the instability is driven by viscous forces in the defending
liquid and resisted by capillary forces at the interface. The stabilising
effect of capillary forces means that sufficiently high wavenumber (i.e. short
wavelength) modes decay; onset occurs as the radius of the interface reaches a
critical value, at which the perimeter of the interface becomes large enough
to accommodate the longest-wavelength unstable mode [5]. In the incompressible
system, when the flow rate is fixed at $\hat{Q}_{b}=1$, Paterson [5] showed
that this critical radius $\hat{s}_{0i}$ satisfies

$\displaystyle\sqrt{\textrm{Ca}\hat{s}_{0i}+\frac{1}{4}}-\frac{1}{2}=2\pi.$
(21)

Hence, $\hat{s}_{0i}\sim\textrm{Ca}^{-1}$, such that onset occurs at smaller
radii for larger Ca or, say, larger $Q$. In the compressible system, the time-
dependent injection rate is $\hat{Q}_{b}=2\hat{s}_{0}\dot{s}_{0}$, where
$\dot{s}_{0}(\hat{t})$ is given by Eq. 7. The critical radius $\hat{s}_{0c}$
in the compressible system must then satisfy

$\displaystyle\sqrt{\frac{2\textrm{Ca}\hat{s}_{0c}}{\mathcal{C}}\left(\frac{\hat{s}_{0c}^{2}-\mathcal{S}^{2}-\hat{t}_{c}}{\ln(\hat{s}_{0c})}\right)+\frac{1}{4}}-\frac{1}{2}=2\pi,$
(22)

where $\hat{t}_{c}$ is the time of onset (i.e.,
$\hat{s}_{0}(\hat{t}_{c})=\hat{s}_{0,c}$). To close Eq. (22), we require
$\hat{s}_{0}(\hat{t})$, which we calculate numerically from the reduced model.
Hence, onset in the compressible system depends on both Ca and $\mathcal{C}$,
with the latter dictating the unsteady injection rate.

Figure 9: Onset radius $\hat{s}_{0c}$ in the compressible system [Eq. (21)],
normalised by the onset radius $\hat{s}_{0i}$ in the incompressible system
[Eq. (22)] at the same Ca, as functions of Ca for $\mathcal{C}$ ranging from
0.01 to 10 (increasing in the direction of the arrows). Red and blue curves,
respectively, indicate $\mathcal{C}\geq 1$ and $\mathcal{C}<1$. The dashed
black line denotes $\hat{s}_{0c}=\hat{s}_{0i}$. Inset: Close-up of the low-Ca
region where compressibility can weakly promote onset
($\hat{s}_{0i}<\hat{s}_{0c}$), according the predictions of linear stability
analysis.

In Figure 9, we compare the linear stability predictions of onset in the
compressible and incompressible systems by plotting
$\hat{s}_{0c}/\hat{s}_{0i}$ as a function of Ca for varying $\mathcal{C}$. We
determine $\hat{s}_{0c}$ numerically, while $\hat{s}_{0i}$ is given
analytically by Eq. (21). We observe a substantial delay in the onset of the
instability, indicated by $\hat{s}_{0c}/\hat{s}_{0i}>1$, when
$\textrm{Ca}\gtrsim 10^{2}$. At the highest Ca and $\mathcal{C}$ shown, onset
is delayed in the compressible system to radii more than 20 times greater than
in the incompressible system at the same Ca. The strong delaying effect
derives from the very low initial injection rates $\hat{Q}_{b}\ll 1$ [see Fig
4(a)] associated with compression-driven displacement. The effect is amplified
with increasing $\mathcal{C}$, which leads to lower initial $\hat{Q}_{b}$, and
with increasing Ca, which corresponds to smaller onset radii in the
incompressible system. Hence, if onset is predicted at small radii
$\hat{s}_{0}\ll 1$ for an incompressible flow, the slow initial injection rate
can significantly delay onset in a compression-driven flow. The inset of Fig.
9 shows a magnified plot of $\hat{s}_{0c}/\hat{s}_{0i}$ at
$\textrm{Ca}\lesssim 10^{3}$. For low Ca, linear stability predicts that
compressibility may act to promote the onset of viscous fingering relative to
an incompressible system, such that $\hat{s}_{0c}/\hat{s}_{0i}<1$. This
promoting effect at low Ca is due to the relatively high injection rates
($\hat{Q}_{b}>1$) at later times or larger $\hat{s}_{0}$, which linear
stability analysis suggests may trigger onset at smaller radii relative to an
incompressible flow. This promoting effect is much weaker than the delaying
effect of high Ca, reducing the radius of onset by a factor of 2. Moreover,
the effect is triggered for large $\hat{s}_{0}\lesssim 1$, and it is therefore
unclear whether the effect on viscous fingering would be visible before
breakout.

###  V.5 Nonlinear finger growth

We now return to fully nonlinear pattern formation in both experiments [Fig.
2] and simulations [Fig. 3]. We first demonstrate that the systematic delay in
the onset of fingering predicted by linear stability analysis due to
compressibility at high Ca is readily observable in both experiments and
simulations. We then examine the impact of this delayed onset on the resulting
fingering pattern. We also examine the agreement between simulations and
experiments, both in terms of the point of onset and the details of the
resulting pattern. To quantify the severity of the fingering pattern, we
consider the isoperimetric ratio

$\displaystyle\mathcal{I}=\frac{L^{2}}{4\pi A},$ (23)

where $L$ and $A$ are, respectively, the length of the interface and the area
it encloses. For a circular interface, $\mathcal{I}=1$. Hence, any deviation
from unity indicates some perturbation away from the axisymmetric base state
assumed in the linear stability analysis [§ V.4].

Figure 10: Isoperimetric ratio $\mathcal{I}$ [Eq. (23)] as a function of the
maximum radial extent of the interface at $\textrm{Ca}=5.21\times 10^{3}$ from
(a) numerical simulations at $\mathcal{C}\in\\{4.5\times 10^{-3},$ 0.009,
0.018, 0.036, 0.071, 0.14, 0.28, 0.57, 1.13, 2.26}, and (b) experiments at
$\mathcal{C}\in\\{0.018$, 0.036, 0.071, 0.14, 0.28} ($Q=2.5$ mL/min,
$V_{g}(0)\in\\{25,50,100,200\\}$ mL). The incompressible solution is shown in
(a) for comparison (dotted black). Dashed curves in (a) are at the values of
$\mathcal{C}$ used for the experiments in (b).

Figure 10 shows the evolution of $\mathcal{I}$ at fixed Ca for a range of
$\mathcal{C}$ for simulations [Figs. 10(a)] and experiments [Figs. 10(b)] as
functions of the maximal radial extent of the inteface $\max(\hat{s})$. These
experiments and simulations correspond to the top row of Figures 2 and 3,
respectively. The isoperimetric ratio $\mathcal{I}$ is initially close to 1,
suggesting that the interface is near-circular, before increasing
monotonically and at a relatively steady rate for the remainder of the
experiment or simulation, corresponding to the growth of viscous fingers. The
departure from $\mathcal{I}\approx 1$ occurs at larger radii [$\max(\hat{s})$]
for larger $\mathcal{C}$ (arrows) in both simulations and experiments,
consistent with a delayed onset. The growth of $\mathcal{I}$ with
$\max(\hat{s})$ is somewhat faster in experiments than in simulations,
consistent with the more frequent occurrence of tip-splitting and side-
branching visible in the experimental images (Fig. 2).

Figure 11: (a) Onset radius $\hat{s}^{*}$ at which $\mathcal{I}$ first exceeds
1.1 (c.f. Fig. 10) and (b) the near-breakout isoperimetric ratio
$\mathcal{I}(\hat{t}_{0.9})$ measured at
$\hat{t}_{0.9}=\hat{t}(\max(\hat{s})=0.9)$ as functions of $\mathcal{C}$ for a
range of Ca. Symbols and colours correspond to $\textrm{Ca}=2.61\times
10^{3}$, $5.21\times 10^{3}$, $1.04\times 10^{4}$, and $2.08\times 10^{4}$
with arrows indicating increasing Ca. Simulation results and experimental
results are plotted with small connected symbols and large scattered symbols,
respectively. The incompressible results are shown for reference (horizontal
dashed lines)

To quantitatively examine how $\mathcal{C}$ impacts the onset of fingering, we
define the radius at onset $\hat{s}^{*}$ in our experiments and simulations as
being $\max(\hat{s})$ at the last recorded instant (video frame or time step)
for which $\mathcal{I}<1.1$. The measured and computed values of $\hat{s}^{*}$
as functions of $\mathcal{C}$ for different Ca are shown in Fig. 11(a) for
experiments (scattered symbols) and simulations (connected symbols).
Qualitatively, both sets of data show the same behaviour: the onset radius
$\hat{s}^{*}$ tends to increase with increasing $\mathcal{C}$ or decreasing
Ca. The latter result is familiar from classical studies of viscous fingering
[5]. Both results are consistent with the linear stability analysis presented
in § V.4, confirming the prediction that increasing $\mathcal{C}$ delays the
onset of fingering. Moreover, varying Ca or $\mathcal{C}$ by a similar amount
(an order of magnitude, say) yields a comparable (and opposite) effect on the
radius of onset. Comparing experiments and simulations directly, there is
reasonable quantitative agreement between the two over the experimental
parameter range studied: the majority of experimental data points lie within
one standard deviation of the corresponding simulation. The gradient of the
experimental data, however, appears shallower than that of the simulations,
which suggests that increasing $\mathcal{C}$ is less effective in delaying
onset in experiments.

The broad range of patterns generated by the fingering instability after onset
can be quantified by considering $\mathcal{I}(\hat{t}_{0.9})$. That is, the
isoperimetric ratio at the near-breakout time
$\hat{t}_{0.9}=\hat{t}(\max(\hat{s})=0.9)$ (see § III.2). We plot
$\mathcal{I}(\hat{t}_{0.9})$ as a function of $\mathcal{C}$ for different
values of Ca in Fig. 11(b). In simulations, we find that
$\mathcal{I}(\hat{t}_{0.9})$ tends to the incompressible case (dashed lines)
as $\mathcal{C}\to 0$, and the severity of the fingering pattern at
$\hat{t}_{0.9}$ decreases with increasing $\mathcal{C}$ or decreasing Ca.
Qualitatively, this behaviour is consistent with experiments, though
experimental measurements are consistently greater, again consistent with the
observed prevalence of tip-splitting and side-branching in experiments (Fig.
2).

Finally, as discussed in § III and § V.2, the systematically greater pressures
observed in experiments than in simulations may be in part due to the
qualitatively and quantitatively different fingering patterns. For instance,
though $\mathcal{I}$ increases more steeply in experiments than in simulations
(Fig. 10), indicating more severe fingering which was previously attributed to
lower pressures, the delayed onset in experiments [Fig. 11(a)] could counter
this effect due to the need to displace more liquid at earlier times.
Furthermore, the greater instances of side-branching and tip-splitting in
experiments may also correspond to greater pressures due to the displacement
of liquid between the primary fingers. The exact correspondence between finger
morphology and driving pressure is beyond the scope of this study, but merits
future investigation.

##  VI Discussion and conclusions

We have studied gas-liquid displacement in a rigid Hele-Shaw cell, driven by
the steady compression of a connected gas reservoir. By considering an
axisymmetric interface, we developed a simple reduced model [Eq. (7)]
analogous to the recent work of Cuttle and MacMinn [29], who studied
compression-driven displacement in a capillary tube. The unsteady injection
rate and the gas pressure in the reduced model are controlled by a single
dimensionless parameter, the compressibility number $\mathcal{C}$, and are
independent of the capillary number Ca. Remarkably, in experiments and
simulations, which are subject to viscous fingering and therefore strongly
non-axisymmetric, we found that the injection rate and gas pressure were still
controlled primarily by $\mathcal{C}$. Variations in Ca had a far more modest
effect on these ‘bulk’ dynamics, despite having a strong influence on the
severity of the fingering instability. We therefore argue that $\mathcal{C}$
is the key control parameter for bulk displacement dynamics, even for
hydrodynamically unstable flows.

The reduced model [Eq. (7)] also revealed two underlying dynamical regimes
that arise from the basic coupling between a viscous displacement flow and the
volumetric compression of a gas. The low- and high-$\mathcal{C}$ regimes
correspond to “on-time” and quasisteady or delayed and burst-like expulsion at
the moment of breakout, when the interface reaches the outlet of the cell. We
rationalised these regimes by following the dynamical-systems approach
employed by Cuttle and MacMinn [29] in studying the corresponding capillary-
tube problem. As in the capillary tube, there exists a critical
compressibility number $\mathcal{C}=1$, where the trivial solutions of the
system vanish for $\hat{s}_{0}\leq 1$; in the tube, the critical point
$\mathcal{C}=1$ dictates the transition between quasi-steady and burst-like
dynamics. However, due to the evolving base state (increasing radius) in the
radial system, we identified a sensitivity to the initial radius that can also
influence whether the dynamics are quasi-steady or burst-like. While the
reduced model neglects the fingering instability, we nonetheless found that
the delayed breakout and over-pressure predicted by the model were
qualitatively recovered in simulations (and experiments; not shown) [Fig. 7],
which again points to the robust role of compression-driven displacement
dynamics in a pattern-forming system.

To understand the impact of compressible displacement dynamics on the onset of
viscous fingering, in § V.4, we performed a linear stability analysis of the
reduced model, which we took as the base state. We found that the growth rate
and the most unstable mode of the perturbations depended strongly on both Ca
and $\mathcal{C}$. Specifically, Ca modulates the relative strengths of the
destabilising viscous and stabilising capillary forces at a given flow rate,
as in the classical system, while $\mathcal{C}$ sets the evolution of the
time-dependent injection rate. Our analysis predicted that, compared to an
incompressible flow at the same nominal injection rate $Q$, compressibility
may either suppress or promote the onset of fingering, depending on Ca. The
promoting effect is relatively weak and would occur at very low Ca, an order
of magnitude lower than the experiments presented in this work, so it is
uncertain whether one could observe its effects in practice. In contrast, the
delaying effect at high Ca is much more pronounced and is indeed readily
observed in both experiments and simulations. In fact, increasing
$\mathcal{C}$ was found to be as effective in delaying onset as decreasing Ca
by a similar factor. The result of this delay is that the severity of the
fingering pattern, as measured by the isoperimetric ratio [Eq. (23)],
decreases substantially as $\mathcal{C}$ increases at fixed Ca.

By conducting extensive experimental and numerical studies in tandem, we were
able to thoroughly compare state-of-the-art simulations with physical data.
The most impressive agreement between the two was found in the volume growth
rate of the bubbles and, to a lesser extent, the evolution of the pressure.
This difference is despite the broad variation in patterns, quantified by the
isoperimetric ratio, observed in experiments and simulations at the same
parameters. Indeed, even our repeat experiments were subject to significant
variability in fingering behavior, and yet were highly reproducible in terms
of injection rate and pressure (Appendix B). These observations speak to the
robust nature of the underlying dynamics of compression-driven displacement
that modulate the growth of the interface, and which can be described to
leading order by the single parameter $\mathcal{C}$.

There are several possible sources for the differences in pattern formation
(e.g., the isoperimetric ratio or the qualitative interface evolution) between
experiments and simulations. One such source is the choice of initial
condition in the simulations [Eq. (13)]. We chose an initial radius $r_{0}$ to
best match with experiments. However, small variations in $r_{0}$, along with
the choice of $\varepsilon$ (initial amplitude of perturbations) and the modes
of perturbations used could have a non-negligible influence on the final shape
of the interface. Further, the experiments are subject small disturbances due
to plate defects, microscopic contaminants, thermal fluctuations, and the
outlet conditions, amongst other culprits, which are not captured by our
model. Despite being small, such disturbances can significantly modify the
fingering pattern, particularly at higher Ca where the interface is far more
distorted. For example, it has been shown that finite perturbations in Hele-
Shaw cells and channels can exert a strong influence on pattern formation [6,
50, 51, 52, 53]. As our simulations and experiments are subject to very
different sources of perturbations, it is unsurprising that they should
produce quantitatively and qualitatively different patterns. These details,
however, do not detract from the key result confirmed by both approaches;
increasing the compressibility number is as effective in delaying onset as
decreasing the capillary number. Compressibility acts passively in two-phase
gas-driven flows and is as natural to the system as viscosity or surface
tension. This is in contrast with previous control strategies [25, 23, 19, 21,
26] that, although effective, are often awkward to implement in practice due
to the restrictions placed on the confining geometry or the choice of fluids.
The compressibility number, meanwhile, is a parameter that can easily be
tuned, without having to alter the system geometry, compliance, or fluid
properties; simply selecting a larger syringe is sufficient. Our results
therefore strongly point to compressibility as a second key parameter in the
assessment and control of viscous fingering in real systems, as discussed in
more detail in a companion study [54]. As a final remark, gas compression and
the associated unsteady flows will continue to be a source of frustration in
many practical and experimental settings, where steady flows are required. For
those wishing to avoid such effects, our study offers a complete framework
that accounts for all relevant parameters.

###### Acknowledgements.

We are grateful to Mr. Clive Baker for technical support. This work was
supported by the European Research Council (ERC) under the European Union’s
Horizon 2020 Programme [Grant No. 805469], by the UK Engineering and Physical
Sciences Research Council (EPSRC) [Grant No. EP/S034587/1], and by the John
Fell Oxford University Press Research Fund [Grant No. 132/012].

##  Appendix A Injection from a compressed gas source

In the main text, we considered injection from a syringe pump, where a fixed
mass of gas is compressed at a constant volumetric rate $Q$. Here, we consider
a second scenario, also commonly used in laboratories: injection from a
compressed gas source via a needle valve. For brevity, we consider a circular
front and neglect thin films, as in the derivation of the reduced model
presented in § II.2.

This injection scenario corresponds to a constant mass flow of air into a
reservoir of fixed volume, such that the number of moles $n(t)$ of gas in the
system increases at a fixed rate $\mathrm{d}n/\mathrm{d}t=\dot{n}$. In
practice, the volumetric rate of gas emitted from the needle valve is
calibrated and set at $Q$ while venting into fixed atmospheric pressure
$p_{\textrm{atm}}$, such that

$Q=\left.\frac{\mathrm{d}V_{g}}{\mathrm{d}t}\right|_{p_{g}=p_{\textrm{atm}}}=\frac{R_{g}T}{p_{\textrm{atm}}}\dot{n},$
(24)

where $R_{g}$ and $T$ are the ideal gas constant and gas temperature,
respectively. However, compression of the gas once it is diverted into the
flow cell means that while $\dot{n}$ may remain constant, the actual flow rate
$Q_{b}$ will no longer be equal to $Q$. From the ideal gas law, we have
$p_{g}(t)V_{g}(t)=[n(0)+\dot{n}t]R_{g}T$, where $V_{g}(t)$ is the total volume
of gas in the reservoir and the cell. Also, $p_{g}(0)V_{g}(0)=n(0)RT$, so the
pressure of the gas is

$p_{g}(t)=\frac{V_{g}(0)(p_{\textrm{atm}}+\gamma\\{\pi/[4s_{0}(0)]+2/b\\}+Qtp_{\textrm{atm}}}{V_{g}(0)+V_{b}(t)-V_{b}(0)},$
(25)

where the initial gas pressure $p_{g}(0)$ is the same as in § II.1. Applying
the non-dimensionalisation of Eq. (2) yields

$\Delta\hat{p}_{g}(\hat{t})=\frac{1}{\mathcal{V}+\hat{s}_{0}^{2}(\hat{t})-\mathcal{S}^{2}}\left\\{\mathcal{V}\left[\mathcal{P}+\frac{1}{\textrm{Ca}}\left(\frac{\pi}{4\mathcal{S}}+2\alpha\right)\right]+\mathcal{P}\hat{t}\right\\}-\mathcal{P}.$
(26)

Finally, taking the limits $\mathcal{V}\gg 1$, and
$\mathcal{P}\gg(\pi/(4\mathcal{S})+2\alpha)/\textrm{Ca}$, we have

$\Delta\hat{p}_{g}(\hat{t})=\frac{2(\hat{t}-\hat{s}_{0}^{2}-\mathcal{S}^{2})}{\mathcal{C}}.$
(27)

Comparing with Eq. (7), we find that in the limit of the reduced model,
compressed gas injection is identical to syringe pump injection. We therefore
expect similar dynamics in both scenarios in practice, as was recently
observed by Peng  _et al._ [20].

##  Appendix B Experimental reproducibility

Figure 12: Experimental reproducibility. (a, b) Interface evolution for
experimental repetitions at $Q=2.5$ mL/min and (left to right)
$V_{\mathrm{res}}(0)\in\\{25,50,100,200\\}$ mL. The experiments in (a)
correspond to the Main Set reported throughout the main document, while those
in (b) are the Repeat Set. (c-e) Non-dimensional volume (c), gauge gas
pressure $\Delta\hat{p}_{g}$ (d), and isoperimetric ratio (e) as functions of
non-dimensional time for the experiments shown in (a; solid lines) and (b;
dotted lines). Arrows indicate increasing $V_{g}(0)$.

Figure 12 shows data from two sets of experiments performed at the same
parameter values ($Q=2.5$ mL/min; $V_{g}(0)$=25, 50, 100, and 200 mL). The
Main Set [Fig. 12(a)] corresponds to the set of data presented throughout the
main text, i.e., in Figs. 2, 4(b, d), 10(b). The Repeat Set [Fig. 12(b)] are
used with the Main Set to calculate the mean values shown in Fig. 11, along
with the experimental error bars.

Comparing the two sets, we see that they both exhibit the same qualitative
response in increasing the initial gas volume $V_{g}(0)$ (arrows), which leads
to the growth of fingers at larger radii. Quantitatively, the dynamical growth
of the bubbles $\hat{V}_{b}(\hat{t})$ [Fig. 12(c)] is remarkably reproducible;
data from the Main Set (solid red curves) and Repeat Set (dotted blue curves)
largely overlap on the scale of the plot. This is despite significant
differences in the evolution of both the gauge gas pressure
$\Delta\hat{p}_{g}$ [Fig. 12(d)] and the viscous fingering pattern, again
quantified by the isoperimetric ratio $\mathcal{I}(\hat{t})$ [Fig. 12(e)].
While the observable onset of fingering (when $\mathcal{I}>1.1$) is in
reasonable agreement between the two sets, the Repeat Set has systematically
greater $\mathcal{I}$ than the Main Set. The Hele-Shaw cell was deconstructed,
cleaned and reassembled between the two sets, which will have introduced small
variations in cell geometry between the sets and is likely the reason for the
systematic offset in $\mathcal{I}(\hat{t})$. The greater variability in both
$\Delta\hat{p}_{g}$ and $\mathcal{I}$ than in $\hat{V}_{b}$ may reflect the
greater sensitivity of pressure to variations in fingering pattern (Fig. 5).

## References

  * Hill [1952] S. Hill, Channeling in packed columns, Chem. Eng. Sci. 1, 247 (1952).
  * Chuoke  _et al._ [1959] R. L. Chuoke, P. Van Meurs, and C. van der Poel, The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media, Trans. 216, 188 (1959).
  * Saffman and Taylor [1958] P. G. Saffman and G. I. Taylor, The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid, Proc. Math. Phys. Eng. Sci. 245, 312 (1958).
  * Homsy [1987] G. M. Homsy, Viscous fingering in porous media, Ann. Rev. Fluid Mech. 19, 271 (1987).
  * Paterson [1981] L. Paterson, Radial fingering in a Hele Shaw cell, J. Fluid Mech. 113, 513 (1981).
  * Couder  _et al._ [1986] Y. Couder, N. Gérard, and M. Rabaud, Narrow fingers in the Saffman-Taylor instability, Phys. Rev. A 34, 5175 (1986).
  * Casademunt [2004] J. Casademunt, Viscous fingering as a paradigm of interfacial pattern formation: Recent results and new challenges, Chaos 14, 809 (2004).
  * Sheng [2010] J. Sheng,  _Modern chemical enhanced oil recovery: theory and practice_ (Gulf Professional Publishing, 2010).
  * Lee  _et al._ [2019] C. Lee, B. Zhao, R. Abouatallah, R. Wang, and A. Bazylak, Compressible-gas invasion into liquid-saturated porous media: Application to polymer-electrolyte-membrane electrolyzers, Phys. Rev. Applied 11, 054029 (2019).
  * Ben-Jacob  _et al._ [2000] E. Ben-Jacob, I. Cohen, and H. Levine, Cooperative self-organization of microorganisms, Adv. Phys. 49, 395 (2000).
  * Matsushita and Fujikawa [1990] M. Matsushita and H. Fujikawa, Diffusion-limited growth in bacterial colony formation, Phys. A: Stat. Mech. Appl. 168, 498 (1990).
  * Jha  _et al._ [2013] B. Jha, L. Cueto-Felgueroso, and R. Juanes, Synergetic fluid mixing from viscous fingering and alternating injection, Phys. Rev. Lett. 111, 144501 (2013).
  * Nicolaides  _et al._ [2015] C. Nicolaides, B. Jha, L. Cueto-Felgueroso, and R. Juanes, Impact of viscous fingering and permeability heterogeneity on fluid mixing in porous media, Water Resour. Res. 51, 2634 (2015).
  * Ben-Jacob  _et al._ [1984] E. Ben-Jacob, N. Goldenfeld, B. Kotliar, and J. Langer, Pattern selection in dendritic solidification, Phys. Rev. Lett. 53, 2110 (1984).
  * Mullins [1963] R. F. Mullins, W. W.and Sekerka, Morphological stability of a particle growing by diffusion or heat flow, J. Appl. Phys. 34, 323 (1963).
  * Park and Homsy [1985] C. Park and G. Homsy, The instability of long fingers in hele–shaw flows, Phys. Fluids 28, 1583 (1985).
  * Dias  _et al._ [2012] E. O. Dias, E. Alvarez-Lacalle, M. S. Carvalho, and J. A. Miranda, Minimization of viscous fluid fingering: A variational scheme for optimal flow rates, Phys. Rev. Lett. 109, 144502 (2012).
  * Morrow  _et al._ [2019] L. C. Morrow, T. J. Moroney, and S. W. McCue, Numerical investigation of controlling interfacial instabilities in non-standard Hele-Shaw configurations, J. Fluid Mech. 877, 1063–1097 (2019).
  * Pihler-Puzović  _et al._ [2012] D. Pihler-Puzović, P. Illien, M. Heil, and A. Juel, Suppression of complex fingerlike patterns at the interface between air and a viscous fluid by elastic membranes, Phys. Rev. Lett. 108, 074502 (2012).
  * Peng  _et al._ [2022] G. G. Peng, C. Cuttle, F. Box, J. H. Guan, A. Juel, C. W. MacMinn, and D. Pihler-Puzović, Trapping and escape of viscous fingers in a soft hele-shaw cell, Phys. Rev. Fluids 7, L062001 (2022).
  * Zheng  _et al._ [2015] Z. Zheng, H. Kim, and H. A. Stone, Controlling viscous fingering using time-dependent strategies, Phys. Rev. Lett. 115, 174501 (2015).
  * Vaquero-Stainer  _et al._ [2019] C. Vaquero-Stainer, M. Heil, A. Juel, and D. Pihler-Puzović, Self-similar and disordered front propagation in a radial hele-shaw channel with time-varying cell depth, Phys. Rev. Fluids 4, 064002 (2019).
  * Al-Housseiny  _et al._ [2012] T. T. Al-Housseiny, P. A. Tsai, and H. A. Stone, Control of interfacial instabilities using flow geometry, Nature Phys. 8, 747 (2012).
  * Bongrand and Tsai [2018] G. Bongrand and P. A. Tsai, Manipulation of viscous fingering in a radially tapered cell geometry, Phys. Rev. E 97, 061101 (2018).
  * Nagatsu  _et al._ [2007] Y. Nagatsu, K. Matsuda, Y. Kato, and Y. Tada, Experimental study on miscible viscous fingering involving viscosity changes induced by variations in chemical species concentrations due to chemical reactions, J. Fluid Mech. 571, 475–493 (2007).
  * Gao  _et al._ [2019] T. Gao, M. Mirzadeh, P. Bai, K. M. Conforti, and M. Z. Bazant, Active control of viscous fingering using electric fields, Nat. Commun. 10, 1 (2019).
  * Sandnes  _et al._ [2012] B. Sandnes, E. G. Flekkøy, and K. J. Måløy, Stick slip displacement of confined granular mixtures: Bubble expansion, Eur. Phys. J. Spec. Top. 204, 19 (2012).
  * Lai  _et al._ [2018] B. Lai, C.-Y. and Rallabandi, A. Perazzo, Z. Zheng, S. E. Smiddy, and H. A. Stone, Foam-driven fracture, Proc. Natl. Acad. Sci. 115, 8082 (2018).
  * Cuttle and MacMinn [2023] C. Cuttle and C. W. MacMinn, Dynamics of compression-driven gas-liquid displacement in a capillary tube, Phys. Rev. Lett. (2023), in press. Preprint: arXiv:2112.12898.
  * Sandnes  _et al._ [2011] B. Sandnes, E. G. Flekkøy, H. A. Knudsen, K. J. Måløy, and H. See, Patterns and flow in frictional fluid dynamics, Nat. Commun. 2, 288 (2011).
  * Park  _et al._ [1984] C.-W. Park, S. Gorell, and G. M. Homsy, Two-phase displacement in hele-shaw cells: experiments on viscously driven instabilities, J. Fluid Mech. 141, 275–287 (1984).
  * Tabeling and Libchaber [1986] P. Tabeling and A. Libchaber, Film draining and the saffman-taylor problem, Phys. Rev. A 33, 794 (1986).
  * Kopf‐Sill and Homsy [1987] A. R. Kopf‐Sill and G. M. Homsy, Narrow fingers in a Hele–Shaw cell, Phys. Fluids 30, 2607 (1987).
  * Tabeling  _et al._ [1987] P. Tabeling, G. Zocchi, and A. Libchaber, An experimental study of the saffman-taylor instability, J. Fluid Mech. 177, 67–82 (1987).
  * Moore  _et al._ [2002] M. G. Moore, A. Juel, J. M. Burgess, W. D. McCormick, and H. L. Swinney, Fluctuations in viscous fingering, Phys. Rev. E 65, 030601 (2002).
  * Ristroph  _et al._ [2006] L. Ristroph, M. Thrasher, M. B. Mineev-Weinstein, and H. L. Swinney, Fjords in viscous fingering: Selection of width and opening angle, Phys. Rev. E 74, 015201 (2006).
  * Maxworthy [1989] T. Maxworthy, Experimental study of interface instability in a hele-shaw cell, Phys. Rev. A 39, 5863 (1989).
  * Martyushev  _et al._ [2009] L. M. Martyushev, A. I. Birzina, M. S. Konovalov, and A. P. Sergeev, Experimental investigation of the onset of instability in a radial hele-shaw cell, Phys. Rev. E 80, 066306 (2009).
  * Chen [1987] J.-D. Chen, Radial viscous fingering patterns in hele-shaw cells, Exps. Fluids 5, 363 (1987).
  * Zhao  _et al._ [2016] B. Zhao, C. W. MacMinn, and R. Juanes, Wettability control on multiphase flow in patterned microfluidics, Proc. Natl. Acad. Sci. 113, 10251 (2016).
  * Li  _et al._ [2009] S. Li, J. S. Lowengrub, J. Fontana, and P. Palffy-Muhoray, Control of viscous fingering patterns in a radial hele-shaw cell, Phys. Rev. Lett. 102, 174501 (2009).
  * Zhao  _et al._ [2020] M. Zhao, P. H. Anjos, J. Lowengrub, and S. Li, Pattern formation of the three-layer saffman-taylor problem in a radial hele-shaw cell, Phys. Rev. Fluids 5, 124005 (2020).
  * Morrow  _et al._ [2023a] L. C. Morrow, N. De Cock, and S. W. McCue, Viscous fingering patterns for hele–shaw flow in a doubly connected geometry driven by a pressure differential or rotation, Phys. Rev. Fluids 8, 014001 (2023a).
  * Peng  _et al._ [2015] G. G. Peng, D. Pihler-Puzović, A. Juel, M. Heil, and J. R. Lister, Displacement flows under elastic membranes. part 2\. analysis of interfacial effects, J. Fluid Mech. 784, 512 (2015).
  * Park and Homsy [1984] C.-W. Park and G. M. Homsy, Two-phase displacement in Hele Shaw cells: theory, J. Fluid Mech. 139, 291 (1984).
  * Morrow  _et al._ [2021] L. C. Morrow, T. J. Moroney, M. C. Dallaston, and S. W. McCue, A review of one-phase Hele-Shaw flows and a level-set method for non-standard configurations, ANZIAM J. 63, 269 (2021).
  * Osher and Sethian [1988] S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79, 12 (1988).
  * Moroney  _et al._ [2017] T. J. Moroney, D. R. Lusmore, S. W. McCue, and D. L. S. McElwain, Extending fields in a level set method by solving a biharmonic equation, J. Comput. Phys. 343, 170 (2017).
  * Gibou  _et al._ [2002] F. Gibou, R. P. Fedkiw, L.-T. Cheng, and M. Kang, A second-order-accurate symmetric discretization of the Poisson equation on irregular domains, J. Comput. Phys. 176, 205 (2002).
  * Zocchi  _et al._ [1987] G. Zocchi, B. E. Shaw, A. Libchaber, and L. P. Kadanoff, Finger narrowing under local perturbations in the Saffman-Taylor problem, Phys. Rev. A 36, 1894 (1987).
  * Rabaud  _et al._ [1988] M. Rabaud, Y. Couder, and N. Gerard, Dynamics and stability of anomalous saffman-taylor fingers, Phys. Rev. A 37, 935 (1988).
  * Thompson  _et al._ [2014] A. B. Thompson, A. Juel, and A. L. Hazel, Multiple finger propagation modes in Hele-Shaw channels of variable depth, J. Fluid Mech. 746, 123 (2014).
  * Zhang  _et al._ [2021] Q. Zhang, A. Amooie, M. Z. Bazant, and I. Bischofberger, Growth morphology and symmetry selection of interfacial instabilities in anisotropic environments, Soft Matter 17, 1202 (2021).
  * Morrow  _et al._ [2023b] L. C. Morrow, C. Cuttle, and C. W. MacMinn, Compression-driven displacement as a passive means of controlling viscous fingering, in review (2023b).