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

TL;DR

This study investigates how gas compressibility influences viscous fingering in a radial Hele-Shaw cell, revealing that higher compressibility delays finger formation and reduces pattern severity, with experiments and simulations providing comprehensive insights.

Contribution

It introduces the compressibility number as a key parameter affecting viscous fingering, expanding understanding beyond traditional incompressible models.

Findings

Increasing compressibility number delays finger onset.

Higher compressibility reduces fingering pattern severity.

Experiments and simulations show consistent effects across parameters.

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 attracted research interest for decades. 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.

Equations

$$\mathbf{v}$$

$${\nabla }^{2}p$$

$$v_n$$

$$p$$

$$f_1=\frac{|v_n\eta /\gamma {|}^{2/3}}{0.76+2.16|v_n\eta /\gamma {|}^{2/3}},\mathrm{and}{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}},$$

$$p_g(t)V_g(t)$$

$$p_g(t)$$

$$\Delta p_g(t)$$

$$\Delta p$$

$$\hat{\mathbf{x}}=\frac{\mathbf{x}}{R_c},\hat{t}=\frac{Q}{\pi R_c^{2}b}t,\hat{p}=\frac{\pi b^3}{6\mu Q}p,\hat{\mathbf{v}}=\frac{2\pi R_{c}b}{Q}\mathbf{v},\hat{V}=\frac{V}{\pi R_c^{2}b},$$

$${\hat{\nabla }}^2\hat{p}$$

$${\hat{v}}_n$$

$$\Delta \hat{p}$$

$$\Delta {\hat{p}}_g$$

$$\text{Ca}=\frac{12{\alpha }^2\mu Q}{2\pi R_{c}b\gamma },\mathcal{C}=\frac{12\mu Q{V}_g(0)}{{\pi }^2{R}_c^2{b}^4{p}_{\text{atm}}},\mathcal{V}=\frac{V_g(0)}{\pi R_c^{2}b},\alpha =\frac{R_c}{b},\mathrm{and}\mathcal{R}=\frac{r_0}{R_c},$$

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

$$\frac{\text{d}{\hat{R}}_0}{\text{d}\hat{t}}$$

$$\Delta {\hat{p}}_g$$

$$\frac{\text{d}{\hat{R}}_0}{\text{d}\hat{t}}$$

$$\frac{\partial \varphi }{\partial \hat{t}}+F\left|\hat{\nabla }\varphi \right|=0,$$

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

$${\hat{V}}_g={\int }_{{\mathbb{R}}^2}(1-f)H(\varphi )\text{d}\hat{V},$$

$$H=\{\begin{array}{ll}0& \text{ if }\varphi <-\delta \\ 1& \text{ if }\varphi >\delta \\ \frac{1}{2}\left[1+\frac{\varphi }{\delta }+\frac{1}{\pi }\sin \left(\frac{\pi \varphi }{\delta }\right)\right]& \text{ otherwise.}\end{array},$$

$$\hat{R}(\theta ,0)=\mathcal{R}\left\{1+\sum _{n=2}^{12}{\epsilon }_n\cos \left[n(\theta -2\pi {\theta }_n)\right]\right\},$$

$$V_b(t)=V_b(0)+Qt-V_{\mathrm{res}}(0)\left[1-\left(\frac{p_{\text{atm}}}{p_{\text{atm}}+\Delta p_g(t)}\right)\right].$$

$$\lambda ({\hat{Q}}_b)=\lambda (\Delta {\hat{p}}_g)-\lambda (\hat{\omega })=\frac{2}{\hat{\omega }\mathcal{C}}\left(\frac{1}{{\hat{Q}}_b}-1+\frac{\mathcal{C}{\hat{Q}}_b}{4{\hat{R}}_0^2}\right),$$

$${\bar{Q}}_{\pm }=2{\hat{R}}_0\left(\frac{{\hat{R}}_0\pm \sqrt{{\hat{R}}_0^2-\mathcal{C}}}{\mathcal{C}}\right),$$

$$\hat{R}(\theta ,\hat{t})$$

$$\hat{p}(\hat{r},\theta ,\hat{t})$$

$$\lambda ({\hat{\gamma }}_n)=\frac{n-1}{2{\hat{R}}_0^2}\left[{\hat{Q}}_b-\frac{n(n+1)}{\text{Ca}{\hat{R}}_0}\right].$$

$$n_{\max }=\sqrt{\frac{1+\text{Ca}{\hat{R}}_0{\hat{Q}}_b}{3}},$$

$$\sqrt{\text{Ca}{\hat{R}}_{0i}+\frac{1}{4}}-\frac{1}{2}=2\pi .$$

$$\sqrt{\frac{2\text{Ca}{\hat{R}}_{0c}}{\mathcal{C}}\left(\frac{{\hat{R}}_{0c}^2-{\mathcal{R}}^2-{\hat{t}}_c}{\ln ({\hat{R}}_{0c})}\right)+\frac{1}{4}}-\frac{1}{2}=2\pi ,$$

Full text

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, 5]. This classical viscous-fingering instability has been extensively studied as an archetype of interfacial pattern formation [6, 7, 8] and for its relevance to enhanced oil recovery [1, 9, 10, 3, 11], for which fingering poses major obstacles. Modern applications include the operation of fuel cells [12, 13], the remediation of groundwater contamination [14, 15], and the subsurface sequestration of CO2 [16, 17] or storage of hydrogen [18]. 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 fluids and opposed by capillary forces at the interface. The capillary number Ca [19], 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 [20, 21], replacing one of the rigid plates with an elastic membrane or slab [22, 23], varying the gap between the plates as the gas is injected [24, 25], using a tapered (rather than uniform) flow cell [26, 27], and the application of external electric fields [28]. Here, we examine a simple but previously overlooked mechanism in classical gas-driven viscous fingering: the compression of the injected gas.

Gas-driven displacement of liquid is a common practical and experimental scenario that is highly susceptible to viscous fingering. In addition to generally 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 [29, 30, 31], which are known to exert a fundamental influence on pattern-forming processes [32]. 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 gas at the inlet [33, 34, 35, 36, 37, 38], performing analysis based on the instantaneous interface velocity [39, 27], or constraining themselves to low injection rates [40]. Alternatively, one may inject a low-viscosity liquid rather than a gas, such as injecting water into oil [41, 42]. 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 the impact of compressibility on displacement flows and viscous fingering in Hele-Shaw cells has not previously been considered in any detail.

Here, we use a combination of mathematical modelling, simulations, and experiments to 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}$ also plays an 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 [31]. 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 [43, 20, 21, 44, 45, 46]. 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.

II Mathematical Model

II.1 Full quasi-2D model

[TABLE]

We consider the system shown schematically in Figure 1, comprising a rigid, circular Hele–Shaw cell of radius $R_{c}$ 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 $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

[TABLE]

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

[TABLE]

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

[TABLE]

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 [48]. Following Peng et al. [47], we model these films via the empirical functions

[TABLE]

which were derived by fitting to simulations of viscous fingering in a rigid Hele-Shaw channel [49]. 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. The function $f_{2}$ captures the modified dynamic boundary condition due to capillary and viscous stresses at the interface. We take the gas to be a fixed mass of ideal gas and the process to be isothermal, meaning that temperature changes due to compression and expansion equilibrate rapidly with the environment. In practice, the compression of the gas will not be perfectly isothermal due to the timescale of thermal diffusion through the glass walls of the syringe. However, our data suggests that an isothermal model is suitable for our experiments, as discussed in Appendix A. Pressure and volume are then related by Boyle’s law,

[TABLE]

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

[TABLE]

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}$ and $V_{b}(0)=\pi r_{0}^{2}b$. While pressure and volume must be related in terms of absolute pressure, we will find it instructive to consider gauge pressures, measured relative to atmospheric pressure. The gauge gas pressure $\Delta p_{g}=p_{g}-p_{\textrm{atm}}$ can be written as

[TABLE]

At the outlet of the Hele–Shaw cell (around the rim), we take the pressure of the liquid to be atmospheric, $p(r=R_{c})=p_{\textrm{atm}}$, such that the gauge liquid pressure $\Delta p=p-p_{\textrm{atm}}$ is

[TABLE]

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

We non-dimensionalize our problem via

[TABLE]

leading to

[TABLE]

where Ca is the capillary number, often referred to as the ‘modified’ capillary number as it incorporates the role of the cell aspect ratio $\alpha$ [e.g., 19]. The parameter $\mathcal{C}$, which we refer to as the compressibility number, is directly analogous to the compressibility number identified in Ref. [31]. Physically, $\mathcal{C}$ may be considered the ratio of viscous and compressive pressure scales or, equivalently, the ratio of the rate of viscous depressurisation (drainage) to that of compressive pressurisation (compression).

For comparison, we also consider the case where the gas is incompressible. In this limit, Eq. (3e) is replaced by the constraint that $\Delta 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{C}$ and $\mathcal{V}$ can be eliminated. This constraint can be enforced by choosing $\Delta p_{g}$ such that

[TABLE]

which ensures that the total rate of outflow is unity.

II.2 Axisymmetric 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=R_{0}(t)$. Equations (3a)-(3e) then reduce to

[TABLE]

We can simplify further by taking $\mathcal{V}\gg 1$ and $2\mathcal{C}^{-1}\gg\textrm{Ca}^{-1}[\pi/(4\mathcal{R})+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 $2\mathcal{C}^{-1}\gg\textrm{Ca}^{-1}[\pi/(4\hat{R}_{0})+2\alpha f_{2}]$ and thus eliminates $f_{2}$ from the model. Additionally ignoring the kinematic impact of thin films (i.e., assuming $f_{1}\approx{}0$), Eqs. (5)-(6) reduce to the ordinary differential equation

[TABLE]

The evolution of the interface then depends on a single dimensionless parameter, the compressibility number $\mathcal{C}$. When the gas is incompressible ($\mathcal{C}\to 0$), Eq. (7) degenerates to $\hat{R}_{0}(\hat{t})=\left(\mathcal{R}^{2}+\hat{t}\right)^{1/2}$ and $\Delta\hat{p}_{g}=\ln(1/\hat{R}_{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 B, we derive a mathematical model for this alternative approach and show that the two methods are identical under the assumptions of the axisymmetric model.

II.3 Numerical scheme

We solve Eqs. (3a)-(3e) using the numerical scheme proposed by Morrow et al. [50], which we briefly summarize here. The scheme is based on the level-set method [51], 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

[TABLE]

where

[TABLE]

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 [52]. To solve Eq. (8), 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|]$. To maintain $\phi$ as a signed distance function such that $|\hat{\nabla}\phi|=1$, we occasionally perform reinitialization [50]. We solve Eq. (3b) for the liquid pressure via a finite difference stencil. Following Gibou et al. [53], 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

[TABLE]

where

[TABLE]

Here, $\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. (1e)]. 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

[TABLE]

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 glass plates of radius $R_{c}=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, 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 $R_{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 C).

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 pressures, 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

[TABLE]

Accounting for $V_{b}(0)$ at this step ensures that the initial value of $V_{b}(t)$ is accurate. In Appendix A we compare Eq. (13) with an adiabatic model to show that the assumption of isothermal compression is justified for our experiments. 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 gas region and the shape of the interface as functions of time. The key steps of the algorithm are illustrated in Fig 2. From each raw image [Fig 2(a)], a reference image of the cell prior to air injection was subtracted. The resulting difference image was then contrast-enhanced using imadjust to facilitate the isolation of the interface [Fig 2(b)]. Pixels along the interface were identified using the edge function with a canny filter [Fig 2(c)]. The resulting binary image was dilated by applying imerode with a circular structuring element of radius one pixel to a negative of the edge-detected image, which ensured connectivity of the interface [Fig 2(d)]. A binary image of the air region was generated by filling the connected interface contour using bwconncomp [Fig 2(e)] and then eroding using imerode with the same structure element to compensate for the earlier dilation step [Fig 2(f)]. 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(R)$ is greater than $0.9R_{c}$, 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 injection rate $Q$ and initial gas volume $V_{g}(0)$ for practical reasons. Moreover, while Ca depends only on $Q$ and not on $V_{g}(0)$, $\mathcal{C}$ is proportional to the product $QV_{g}(0)$. To vary $\mathcal{C}\propto QV_{g}(0)$ 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)$. 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{V}$ are given in the caption of each figure, where appropriate. The remaining parameters are fixed at $\alpha=250$ and $\mathcal{R}=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 injection rate and injection pressure, that arise from compression-driven displacement in a radial Hele-Shaw cell. We show that the axisymmetric model derived in § II.2 qualitatively captures the variation in bulk displacement dynamics with increasing $\mathcal{C}$. Comparison between experiments, simulations and the axisymmetric 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 injection rate in the axisymmetric 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 may 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

We begin by briefly examining the qualitative impact of compressibility on viscous fingering. In Figs. 3 and 4, we show the evolution of the fingering pattern across all of our experiments and simulations, respectively. For the experiments [Fig. 3], 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 fixed along each north-west/south-east diagonal and increasing from left to right in each row. The behaviour for increasing Ca (top to bottom) is as expected from classical work on viscous fingering [41]: 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.

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 [32], 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 [31]. 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 4, 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. 3. 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.

V.2 Bulk displacement dynamics

The dynamics of the axisymmetric model [Eq. (7)] are illustrated in Fig. 5, which shows the evolution of the nondimensional injection rate $\hat{Q}_{b}$ [Fig. 5(a)] and the gauge gas pressure $\Delta\hat{p}_{g}$ [Fig. 5(c)] for different compressibility numbers $\mathcal{C}$. For reference, the incompressible solution for a circular interface is also shown (dashed black curves), for which $\hat{Q}_{b}=1$ and $\Delta\hat{p}_{g}$ monotonically decreases as liquid drains and viscous resistance decreases. In the axisymmetric 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{R}_{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 [31], as discussed further in § V.3.

We show in Fig. 5(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 axisymmetric 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 axisymmetric 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 axisymmetric 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. 3 and 4, 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).

As shown in Fig. 6, the evolution of $\hat{Q}_{b}$ [Fig. 6(a)] and $\Delta\hat{p}_{g}$ [Fig. 6(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. 3 and fifth column of Fig. 4). 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. 5(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. 3 and 4] may also influence the pressure evolution, as illustrated by the results of Fig. 6(b) (see § V.5).

V.3 Dynamical systems framework

The axisymmetric model exhibits two dynamical regimes, characterised by the change from non-monotonic to monotonically increasing $\hat{Q}_{b}$ around $\mathcal{C}\approx 1$ [Fig. 5$(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 [31]; we take the same approach here to describe the dynamics of axisymmetric displacement in a Hele-Shaw cell, as embodied by the axisymmetric model [Eq. (7)].

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

[TABLE]

where $\lambda(\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 axisymmetric model admits two trivial solutions, for which the driving compressive force and the opposing viscous resistance decrease at the same relative rate [i.e., $\lambda(\hat{Q}_{b})=0$]. These are

[TABLE]

which satisfy $\lambda(\bar{Q}_{\pm})=0$.

The influence of these trivial solutions on compressible displacement dynamics is illustrated by the phase-space plot in Fig. 7, where $\bar{Q}_{\pm}(\hat{R}_{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 axisymmetric model is plotted as a faint blue curve, displaying the characteristic nonmonotonic variation in $\hat{Q}_{b}$ with $\hat{R}_{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{R}_{0}$-dependence, the trivial solutions only exist for $\hat{R}_{0}\geq\sqrt{\mathcal{C}}$. At $\hat{R}_{0}=\sqrt{\mathcal{C}}$, the two branches $\bar{Q}_{\pm}$ meet and annihilate at a saddle-node bifurcation. For $\hat{R}_{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{R}_{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{R}$. Computing the basin of attraction in terms of $\mathcal{R}$ for a given $\mathcal{C}$ is only possible numerically, and is beyond the scope of this study. A brief exploration suggests that, for $\mathcal{R}=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 [31].

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

In Figures 8(a, b), we plot the breakout time and the breakout pressure, respectively. We compute these values numerically from the axisymmetric model (solid black curves) for $\mathcal{R}=0.025$. We find that, for $\mathcal{C}\lesssim 1$, the breakout time $\hat{t}_{f}$ is almost exactly 1 (to within numerical resolution) while the breakout pressure $\Delta\hat{p}_{g}(\hat{t}_{f})$ is almost exactly zero (to within numerical resolution). This corresponds to a scenario where breakout occurs at exactly the same time as for an incompressible flow driven at the nominal injection 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. 8 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 axisymmetric 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

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

[TABLE]

where $\varepsilon\ll 1$, $n\geq 2$, and $\hat{R}_{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 [6], we use the $\mathcal{O}(\varepsilon)$ problem to derive an evolution equation for the relative growth rate of the $n$th mode of perturbation,

[TABLE]

The most unstable mode of perturbation is then

[TABLE]

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

Figures 9(a, c) show $\lambda(\hat{\gamma}_{\max})$ and $n_{\max}$, respectively, as functions of $\hat{R}_{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{R}_{0})$ into Eqs. (18) and (19). Increasing Ca enhances the growth rate of the instability [$\lambda(\hat{\gamma}_{n})$ increases] as well as the most unstable mode of perturbation $n_{\max}$ for all $\hat{R}_{0}$. These observations are consistent with our experimental and numerical results (Figs. 3 and 4) where, over the values of $\mathcal{C}$ considered, increasing Ca results in more prominent branching and tip splitting behaviour.

Fixing Ca and varying $\mathcal{C}$, as in Figs. 9(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 $\lambda(\hat{\gamma}_{n})$ and $n_{\max}$ cross the incompressible solution (dashed line; $\hat{Q}_{b}=1$) at some value of $\hat{R}_{0}$.

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 [6]. In the incompressible system, where the injection rate is fixed at $\hat{Q}_{b}=1$, Paterson [6] showed that this critical radius $\hat{R}_{0i}$ satisfies

[TABLE]

Hence, $\hat{R}_{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{R}_{0}\dot{R}_{0}$, where $\dot{R}_{0}(\hat{t})$ is given by Eq. 7. The critical radius $\hat{R}_{0c}$ in the compressible system must then satisfy

[TABLE]

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

In Figure 10, we compare the linear stability predictions of onset in the compressible and incompressible systems by plotting $\hat{R}_{0c}/\hat{R}_{0i}$ as a function of Ca for varying $\mathcal{C}$. We determine $\hat{R}_{0c}$ numerically, while $\hat{R}_{0i}$ is given analytically by Eq. (20). We observe a substantial delay in the onset of the instability, indicated by $\hat{R}_{0c}/\hat{R}_{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 5(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{R}_{0}\ll 1$ for an incompressible flow, then the slow initial injection rate can significantly delay onset in a compression-driven flow. The inset of Fig. 10 shows a magnified plot of $\hat{R}_{0c}/\hat{R}_{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{R}_{0c}/\hat{R}_{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{R}_{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 less than a factor of 2. Moreover, the effect is triggered for large $\hat{R}_{0}\lessapprox 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. 3] and simulations [Fig. 4]. 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

[TABLE]

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 11 shows the evolution of $\mathcal{I}$ at fixed Ca for a range of $\mathcal{C}$ for simulations [Figs. 11(a)] and experiments [Figs. 11(b)] as functions of the maximal radial extent of the interface $\max(\hat{R})$. These experiments and simulations correspond to the top row of Figures 3 and 4, 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{R})$] for larger $\mathcal{C}$ (arrows) in both simulations and experiments, consistent with a delayed onset. The growth of $\mathcal{I}$ with $\max(\hat{R})$ 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. 3).

To quantitatively examine how $\mathcal{C}$ impacts the onset of fingering, we define the radius at onset $\hat{R}^{*}$ in our experiments and simulations as being $\max(\hat{R})$ at the last recorded instant (video frame or time step) for which $\mathcal{I}<1.1$. The measured and computed values of $\hat{R}^{*}$ as functions of $\mathcal{C}$ for different Ca are shown in Fig. 12(a) for experiments (scattered symbols) and simulations (connected symbols). Qualitatively, both sets of data show the same behaviour: the onset radius $\hat{R}^{*}$ tends to increase with increasing $\mathcal{C}$ or decreasing Ca. The latter result is familiar from classical studies of viscous fingering [6]. 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{R})=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. 12(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 systematically greater, again consistent with the observed prevalence of tip-splitting and side-branching in experiments (Fig. 3).

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. 11), indicating more severe fingering which was previously attributed to lower pressures, the delayed onset in experiments [Fig. 12(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.

Finally, we return to the bulk displacement dynamics to shed some light on the increasingly delayed onset and mitigated fingering with $\mathcal{C}$. The near-breakout volume of air in the cell, $\hat{V}_{b}(\hat{t}_{0.9})$, which is approximately equal to the near-breakout volume of liquid expelled since $\hat{V}_{b}(0)\ll{}\hat{V}_{b}(\hat{t}_{0.9})$, is close to the incompressible value and essentially independent of $\mathcal{C}$ for $\mathcal{C}\lesssim 0.1$ [Fig. 13(a)]. Recall that the same is true of the breakout time $\hat{t}_{f}$ [Fig. 8(a)]. Thus, the average injection rate up to near-breakout (i.e., the average value of $\hat{Q}_{b}$ during the interval $0\leq\hat{t}\leq\hat{t}_{0.9}$), which is given by $\hat{V}_{b}(\hat{t}_{0.9})/\hat{t}_{0.9}$, is approximately unity for $\mathcal{C}\lesssim 0.1$ [Fig. 13(b)]. That is, the average injection rate is approximately equal to the nominal injection rate for $\mathcal{C}\lesssim 0.1$, despite the non-trivial time evolution of the instantaneous injection rate $\hat{Q}_{b}(t)$. Nonetheless, these smaller values of $\mathcal{C}$ lead to noticeably delayed onset and reduced intensity of fingering [Fig. 12]. Both $\hat{V}_{b}(\hat{t}_{0.9})$ and $\mathcal{I}(\hat{t}_{0.9})$ increase substantially as $\mathcal{C}$ increases further ($\mathcal{C}\gtrsim{}0.1$) [Figs. 12b and 13b], but $\hat{V}_{b}(\hat{t}_{0.9})/\hat{t}_{0.9}$ decreases substantially [Fig. 13]. Thus, weak compressibility ($\mathcal{C}\lesssim{}0.1$) delays onset and mitigates fingering by introducing a time-varying $\hat{Q}_{b}(\hat{t})$ while roughly preserving the average injection rate, whereas stronger compressibility ($\mathcal{C}\gtrsim{}0.1$) further delays onset, mitigates fingering, and increases the volume of liquid expelled by reducing the average injection rate.

Note that we observe systematically lower values of $\hat{V}_{b}(\hat{t}_{0.9})$ in experiments than in simulations [Fig. 13a], which is consistent with the systematically larger values of $\mathcal{I}(\hat{t}_{0.9})$ [Figs. 12b] and further suggests that the experiments are subject to more severe fingering than the simulations at the same values of $\mathrm{Ca}$ and $\mathcal{C}$.

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 axisymmetric model [Eq. (7)] analogous to the recent work of Cuttle and MacMinn [31], who studied compression-driven displacement in a capillary tube. The unsteady injection rate and the gas pressure in the axisymmetric 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 axisymmetric 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 [31] in studying the corresponding capillary-tube problem. As in the capillary tube, there exists a critical compressibility number $\mathcal{C}=1$, which dictates the transition between quasi-steady and burst-like dynamics. In our axisymmetric model, $\mathcal{C}$ plays a directly analogous role, with the trivial solutions of the system vanishing for $\hat{R}_{0}\leq 1$ at $\mathcal{C}=1$. However, due to the evolving base state (increasing radius) of the axisymmetric model, we identified an additional sensitivity to the initial radius $\mathcal{R}$ that can influence whether the dynamics are quasi-steady or burst-like. While the axisymmetric 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. 8], which again points to the robust role of compression-driven displacement dynamics in a pattern-forming system.

To understand the impact of compression-driven displacement dynamics on the onset of viscous fingering, in § V.4, we performed a linear stability analysis of the axisymmetric 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. (22)], decreases substantially as $\mathcal{C}$ increases at fixed Ca. This mitigated finger growth can be largely attributed to the delayed onset predicted by linear stability analysis.

In the context of viscous fingering with incompressible fluids, numerous studies have considered the impact of imposing a time-varying injection rate $Q_{b}(t)$, typically with the goal of identifying the $Q_{b}(t)$ that minimises or otherwise controls the number of fingers that develop. For example, one pair of studies considered strategies to suppress fingering by varying $Q_{b}(t)$ while keeping the average injection rate constant (i.e., while still injecting a given total volume in the same time total time) [54, 20]. Dias et al. [54] found that a piecewise-constant $Q_{b}(t)$ with a small initial rate followed by a larger subsequent rate was effective in suppressing onset, whereas Dias et al. [20] found that the optimal form of $Q_{b}(t)$ was linearly increasing in time. Although compressibility leads to a natural and passive (rather than actively controlled) variation in the injection rate, our results share several features with these previous works. Specifically, the time-varying rates $Q_{b}(t)$ observed here for $\mathcal{C}\lesssim 0.1$ [Fig. 5a,b] mimic a small-to-large variation that preserves the average rate, as suggested by Dias et al. [54] [Fig. 13(b)]. In addition, the rates observed here for $\mathcal{C}\gtrsim 0.1$ mimic the monotonically increasing rates suggested by Dias et al. [20]. Thus, although we have not specifically investigated optimisation here, our results suggest that compression-driven displacement at $\mathcal{C}\approx 0.1$ would passively achieve the strongest delay in onset while preserving the nominal injection rate on average, and would therefore be optimal (in the sense of [54]) at a given 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 behaviour, and yet were highly reproducible in terms of injection rate and pressure (Appendix C). 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. (12)]. 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 [7, 55, 56, 57, 58]. 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 [59, 26, 22, 24, 28] 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 [60]. 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.

Data availability — The supporting data for this study are openly available on Zenodo [61].

Published version — This article is published as C. Cuttle, L. C. Morrow, and C. W. MacMinn. Compression-driven viscous fingering in a radial Hele-Shaw cell. Physical Review Fluids, 8:113904, 2023.

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 Isothermal and adiabatic models of gas compression

In our experiments, we assumed that gas compression was isothermal, corresponding to ‘rapid’ equilibration with the environment due to heat diffusion through the walls of the syringe. In practice, the syringe walls are borosilicate glass, which is a poor thermal conductor, and our assumptions therefore require some justification. For instance, one may calculate the typical timescale of thermal diffusion $\tau_{D}=W^{2}/\alpha_{g}\approx 7$ s, based on the thickness of the syringe walls $W=2.2$ mm and the thermal diffusivity of borosilicate glass $\alpha_{g}\approx 0.6$ mm2/s. Alternatively, if we consider thermal diffusion in the air, which has a diffusivity of $18$ mm2/s and a lengthscale on the order of the syringe barrel inner radius $16.3$ mm, we compute a similar timescale of around 14 s. We can then compute the Fourier number $\mathrm{Fo}=t_{f}/\tau_{D}$, which compares the total time $t_{f}$ of an experiment to the timescale of thermal diffusion. Hence, when $\mathrm{Fo}\gg 1$, we expect any heat generated to dissipate rapidly to the environment and maintain a near-constant temperature during compression (i.e., an isothermal process). We find that, for $Q=[1.25,2.5,5,10]$ mL/min, $\mathrm{Fo}\approx[40,20,10,5]$, with some scatter around these mean values due to experimental variability and compressibility effects. Hence, we expect compression to be approximately isothermal for all but the fastest experiments (largest $Q$), where $\mathrm{Fo}\gtrsim 1$ and compression may result in non-negligible changes in temperature.

Rather than assume the gas compression is isothermal, we may have instead assumed that the compression is adiabatic, such that the system is perfectly insulated and no heat is conducted to the environment. (In reality, the process will lie somewhere between these two extremes, but a mixed model is unnecessarily complicated and would cloud the key results of our study.) To compare the two assumptions, we can compute the gas volume from the recorded gas gauge pressure (as discussed in § III.2) for an adiabatic process, for which $p_{g}V_{g}^{\eta}$ is a constant. The adiabatic index $\eta=7/5$ for diatomic gases, of which air is predominantly composed. Equation (13) then becomes

[TABLE]

Figure 14 compares isothermal (13) and adiabatic (23) calculations of the normalised bubble volume $\hat{V}_{b}$ as functions of normalised time $\hat{t}$ for the smallest nominal injection rate $Q$ and reservoir volume $V_{\mathrm{res}}(0)$ [Fig. 23(a)] and the largest $Q$ and $V_{\mathrm{res}}(0)$ [Fig. 23(b)], corresponding to the smallest and largest compressibility numbers, respectively. On the one hand, when $\mathcal{C}\propto QV_{\mathrm{res}}(0)$ is smaller, the difference between adiabatic and isothermal models is minimal [Fig. 23(a)] because the pressure $\Delta p_{g}$ required to drive the flow is smaller and the volume of the bubble $V_{b}$ is more comparable to that of the reservoir $V_{\mathrm{res}}$. Thus, relative changes in the total volume of the gas are smaller and have less impact on the volume of gas in the cell. On the other hand, when $\mathcal{C}$ is larger, the two models predict significantly different $\hat{V}_{b}$ for the same $\Delta p_{g}$ [Fig. 14(b)]. Because $\mathrm{Fo}\gtrsim 1$ for the largest $Q$, it is unclear which of the models is most suitable. However, we can compare the data to an independent upper bound on $V_{b}$, which we calculate from images of the expanding fingering pattern, multiplying the projected area of the pattern $A(t)$ by the depth of the cell $b$. The projected volume $A(t)b$ is then the maximum volume that the air could occupy in the cell at time $t$ in the absence of residual films. Comparing predictions of $\hat{V}_{b}$ from the isothermal and adiabatic models against the normalised projected volume $\hat{A}=Ab/(\pi R_{c}^{2}b)$ (black dashed line), we see in Fig. 14(b) that only the isothermal model stays consistently below this upper bound, while the adiabatic model prediction far exceeds it. Hence, the isothermal model is the better choice for our experiments, even for modest $\mathrm{Fo}$.

Appendix B 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 axisymmetric 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

[TABLE]

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 injection 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

[TABLE]

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

[TABLE]

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

[TABLE]

Comparing with Eq. (7), we find that in the limit of the axisymmetric 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. [23].

Appendix C Experimental reproducibility

Figure 15 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. 15(a)] corresponds to the set of data presented throughout the main text, i.e., in Figs. 3, 5(b, d), 11(b). The Repeat Set [Fig. 15(b)] are used with the Main Set to calculate the mean values shown in Fig. 12, along with the experimental error bars.

Comparing the two sets, we see that they both exhibit the same qualitative response to 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. 15(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. 15(d)] and the viscous fingering pattern, again quantified by the isoperimetric ratio $\mathcal{I}(\hat{t})$ [Fig. 15(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. 6).