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

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.


I. INTRODUCTION
When a fluid is displaced from a Hele-Shaw cell or porous medium by the injection of a less viscous fluid, the fluidfluid 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,3,[9][10][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 CO 2 [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 [27,39], 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, FIG.1: Schematic diagram of the fingering problem considered here.A reservoir of gas (white) is compressed at a constant volumetric rate Q in order to displace liquid (blue) from a Hele-Shaw cell comprising the uniform gap b between two plates.
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 C that compares the rate of viscous depressurisation to the rate of compressive pressurisation.We find that 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 C is increased.We further show that increasing 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 [20,21,[43][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 A. Full quasi-2D model
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 µ.Gas is injected through a hole in the centre of the cell by compressing an attached gas reservoir of initial volume V 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 Ω 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 v to its pressure p via Further, the incompressibility of the liquid phase requires that ∇ • v = 0 and, hence, that Following Peng et al. [47], we impose the kinematic and dynamic conditions at the moving interface ∂Ω(t): where v n is the normal velocity of the interface, n is the unit normal to the interface, γ is the interfacial tension, and κ 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 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, where V g (t) = V b (t) + V 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 res the volume of the gas reservoir.Due to the imposed steady compression rate Q, the volume of gas in the reservoir is V res (t) = V res (0) − Qt, such that the pressure of the gas is 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 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 C) such that the initial in-plane curvature of the interface is 1/r 0 and V b (0) = πr 2 0 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 ∆p g = p g − p atm can be written as 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 atm , such that the gauge liquid pressure ∆p = p − p atm is The flow is driven by the gauge liquid pressure difference between the liquid-gas interface and the outlet.We non-dimensionalize our problem via where The thin-film factors f 1 and f 2 are still given by Eq. (1e), which is dimensionless.The dimensionless control parameters are where Ca is the capillary number, often referred to as the 'modified' capillary number as it incorporates the role of the cell aspect ratio α [e.g., 19].The parameter C, which we refer to as the compressibility number, is directly analogous to the compressibility number identified in Ref. [31].Physically, 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 ∆p g must be chosen so that the actual injection rate Qb ≡ d Vb /d t = 1 at all times and the parameters C and V can be eliminated.This constraint can be enforced by choosing ∆p g such that which ensures that the total rate of outflow is unity.

B. 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 We can simplify further by taking V ≫ 1 and 2C −1 ≫ Ca −1 [π/(4R) + 2α]; 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 2C −1 ≫ Ca −1 [π/(4 R0 ) + 2αf 2 ] and thus eliminates f 2 from the model.Additionally ignoring the kinematic impact of thin films (i.e., assuming f 1 ≈ 0), Eqs. ( 5)-( 6) reduce to the ordinary differential equation The evolution of the interface then depends on a single dimensionless parameter, the compressibility number C. When the gas is incompressible (C → 0), Eq. ( 7) degenerates to R0 ( t) = R 2 + t 1/2 and ∆p g = ln(1/ R0 ), 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.

C. 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 ϕ such that ϕ < 0 in the gas region and ϕ > 0 in the liquid region.We evolve ϕ via the level-set equation where and n = ∇ϕ/| ∇ϕ| 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 x ∈ R 2 \Ω( 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-variationdiminishing Runge-Kutta with ∆t = ∆x/[4 max |F |].To maintain ϕ as a signed distance function such that | ∇ϕ| = 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 κ = ∇ • n.Further, the volume of the gas region is computed with where Here, δ = 1.5∆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 where θ n and ε n are selected at random from uniform distributions on the intervals (0, 1) and (5 × 10 −4 , 10 −3 ), respectively.Simulations are performed on the computational domain 0 ≤ r ≤ 1 and 0 ≤ θ < 2π using 1000 × 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 t = tf .

A. 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 ± 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 B and § V E).
The cell is initially filled with 1,000 cSt silicone oil (dynamic viscosity µ = 0.97 Pa s and surface tension γ = 21 mN m −1 at laboratory temperature 22 ± 1 • 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 ± 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 res (0) ∈ {25, 50, 100, 200} mL, with a relative error of 2-4%.To achieve V res (0)=200 mL, we connect an additional sealed acrylic box of internal volume 120 ± 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 ± 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 res (0).
To conduct an experiment, the air was then compressed using a syringe pump (AL-4000; WPI) at a fixed volumetric rate Q ∈ {1.25, 2.50, 5.00, 10.0} mL min −1 .The gauge pressure ∆p g of the air relative to atmospheric pressure p 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 µ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 E and Appendix C).

B. 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 atm + ∆p g (t)]V g (t) = p atm V res (0).Taking V g (0) = V res (0) and p g (0) = p 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 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 = dV b /dt (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 χ 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 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 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), C is proportional to the product QV g (0).To vary C ∝ QV g (0) while keeping Ca ∝ Q fixed, we fix Q and vary the initial gas volume within the range V g (0) ∈{3.125, 6.25, 12.5, 25, 50, 100, 200, 400, 800, 1600} mL, equivalent to V ∈{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 C fixed, we vary Q while fixing the product QV g (0).We restrict our simulations to the experimental values of Q, corresponding to Ca ∈ {2.61 × 10 3 , 5.21 × 10 3 , 1.04 × 10 4 , 2.08 × 10 4 }.Values of V are given in the caption of each figure, where appropriate.The remaining parameters are fixed at α = 250 and R = 0.025, with dimensional values listed in § III A.

V. RESULTS
Section V is organised as follows.We begin in § V A by describing the key observations of compression-driven viscous fingering.In § V B, 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 B qualitatively captures the variation in bulk displacement dynamics with increasing C. Comparison between experiments, simulations and the axisymmetric model demonstrates that the bulk displacement dynamics of the full system are also controlled primarily by C and are remarkably insensitive to both variations in Ca and the presence of viscous fingering.In § V C, we examine the underlying dynamical-systems structure that dictates the unsteady injection rate in the axisymmetric model.The dynamical regimes described in In § V D, 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 E, we examine the growth of viscous fingers in both experiments and simulations.We observe a strong delay in the onset with increasing 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.

A. 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 nearbreakout time at which the experiment is concluded ( § III B).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 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 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 res (0) does not change Ca, so its impact on viscous fingering is not considered in classical studies.Rather, changing V res (0) changes how 'compressible' the system is [32], as measured by the value of 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 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 C are varied.We analyse and compare these results quantitatively in § V E.

B. 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 Qb [Fig.5(a)] and the gauge gas pressure ∆p g [Fig.5(c)] for different compressibility numbers C. For reference, the incompressible solution for a circular interface is also shown (dashed black curves), for which Qb = 1 and ∆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 Qb and ∆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 ( Qb > 1), at which point the gas begins to expand and ∆p g decreases.For C ≪ 1, the compressible dynamics differ only weakly from the incompressible dynamics, with Qb rapidly reaching and then exceeding the nominal flux, before relaxing back toward the incompressible solution Qb = 1.As C approaches unity, the maximum in Qb increases and occurs at later times.For C > 1, the maximum in Qb vanishes, and the injection rate instead increases monotonically, diverging as the interface escapes the cell at the moment of breakout (when R0 = 1).The qualitative change in dynamics around C = 1 is consistent with the distinct dynamical regimes observed for compression-driven displacement in a capillary tube [31], as discussed further in § V C.
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 Qb and ∆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 Qb occurs at a slightly lower value of C in the simulations and experiments than in the axisymmetric model.Otherwise, the evolution of Q b and ∆p g and the variation with 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 E).
As shown in Fig. 6, the evolution of Qb [Fig.6(a)] and ∆p g [Fig.6(b)] at fixed 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 ∆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 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 FIG. 6: Capillary number Ca has a modest impact on the evolution of injection rate Qb and gauge gas pressure ∆p g .Here, we plot (a) Qb and (b) ∆p g for C = 0.14 from numerical simulations of viscous fingering at Ca ∈ {2.61 × 10 3 , 5.21 × 10 3 , 1.04 × 10 4 , 2.08 × 10 4 } (V ∈ {1.72, 3.44, 6.87, 13.7}; column 5 in Fig. 4) and from the axisymmetric model (dashed black).

C. Dynamical systems framework
The axisymmetric model exhibits two dynamical regimes, characterised by the change from non-monotonic to monotonically increasing Qb around C ≈ 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 Qb = 2 R0 (d R0 /d t).From the axisymmetric model, Qb = ∆p g /ω, where we have introduced the resistance ω = ln 1/ R0 by analogy with Ohm's law.We can then write the axisymmetric model as where λ(x) = ẋ/x is the relative rate of change of the variable x, with ẋ = dx/d 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., λ( Qb ) = 0].These are which satisfy λ( Q± ) = 0.The influence of these trivial solutions on compressible displacement dynamics is illustrated by the phase-space plot in Fig. 7, where Q± ( R0 ) are plotted as dashed and solid black curves, respectively, for C = 0.25.The trivial solutions separate regions where Qb is positive (orange) or negative (purple).The local change of sign in Qb determines the stability of each trivial solution.Specifically, Q+ is a repeller and 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 Qb with R0 and showing that, rather than tending back toward the incompressible solution ( Qb = 1), the system is drawn onto the attractive solution Q− .
Due to their R0 -dependence, the trivial solutions only exist for R0 ≥ √ C. At R0 = √ C, the two branches Q± meet and annihilate at a saddle-node bifurcation.For R0 < √ C there are no trivial solutions: Qb > 0 for all Qb , so that the interface accelerates monotonically.Because R0 must be less than unity, C > 1 implies that there is no attractive solution and the interface accelerates monotonically toward breakout.Even for C < 1, however, the dynamics may fail to converge onto Q− , depending on the initial radius R. Computing the basin of attraction in terms of R for a given C is only possible numerically, and is beyond the scope of this study.A brief exploration suggests that, for space for which Qb is positive or negative are shaded orange or purple, respectively (see colorbar).When λ(∆p g ) > λ(ω), the interface accelerates, Qb > 0. This occurs in three distinct subregions (see inset): when Qb < 1 because the gas is compressing while resistance decreases (∆ ṗg > 0 and ω < 0); when 1 < Qb < Q− because the gas is decompressing slowly (∆ ṗg < 0 is small); and when Qb > Q+ because the gas is over-compressed (∆p g is large).When λ(∆p g ) < λ(ω), the interface decelerates, Qb < 0. This occurs only in the region Q− < Qb < Q+ , where the gas decompresses quickly (∆ ṗg < 0 and |∆ ṗg | is large).R = 0.01, the system fails to converge onto Q− when C ≥ 0.92.Therefore, 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 C = 1 is a precise indicator of the flow regime [31].
The stable trajectories that do converge onto Q− terminate with a breakout flux of Qb ( R0 = 1) ≡ 2 1 − √ 1 − C /C because the resistance ω vanishes at the moment of breakout, driving the system exactly onto Q− .Similarly, for unstable trajectories that fail to converge onto Q− , the vanishing resistance drives divergent λ( Qb ), and hence divergent Qb , 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 R = 0.025.We find that, for C ≲ 1, the breakout time tf is almost exactly 1 (to within numerical resolution) while the breakout pressure ∆p g ( tf ) 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 Q− .For C ≳ 1, in contrast, breakout is delayed ( tf > 1) which means that the air is still compressed at breakout and ∆p g ( tf ) > 0. Hence, the driving pressure remains finite as the opposing resistance vanishes and Qb diverges.
For comparison, in Fig. 8 we also show the breakout time and the breakout pressure from simulations with fingering over a range of C at different values of Ca (colored curves and symbols).We observe qualitatively similar behavior, with tf and ∆p g ( tf ) varying slowly for C ≲ 0.1 before increasing sharply around C ≈ 0.1.We also observe nonmonotonic variation in ∆p g ( tf ) for the highest Ca studied, and we would expect to see the same for the lowest Ca were the range of 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.

D. 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 R(θ, t) = R0 ( t) + εγ n ( t) cos nθ + O(ε 2 ), (16) p(r, θ, t) = p0 (r, t) + ε Ân (r, t) cos nθ + O(ε 2 ), (17) where ε ≪ 1, n ≥ 2, and R0 and p0 are the unperturbed circular solution (i.e., the base state).Here, γn and Ân denote, respectively, the amplitudes of the nth mode of perturbation to the radius and pressure.Following Paterson [6], we use the O(ε) problem to derive an evolution equation for the relative growth rate of the nth mode of perturbation, The most unstable mode of perturbation is then which comes about by solving ∂λ(γ n )/∂n = 0. Thus, the value of n max depends on the evolution of Qb ( R0 ) for an unperturbed circular interface.This base state is precisely the solution to the axisymmetric model, where Qb = 2 R0 Ṙ0 .
By combining the axisymmetric model and linear stability analysis, we may therefore understand how compressiondriven displacement modifies the onset of viscous fingering.Figures 9(a, c) show λ(γ max ) and n max , respectively, as functions of R0 for fixed C at different values of Ca.These curves are calculated by solving Eq. ( 7) numerically and substituting Qb ( R0 ) into Eqs.( 18) and (19).Increasing Ca enhances the growth rate of the instability [λ(γ n ) increases] as well as the most unstable mode of perturbation n max for all R0 .These observations are consistent with our experimental and numerical results (Figs. 3 and 4) where, over the values of C considered, increasing Ca results in more prominent branching and tip splitting behaviour.Fixing Ca and varying C, as in Figs.9(b, d), we find that increasing 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 B, compressiondriven 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 C, both λ(γ n ) and n max cross the incompressible solution (dashed line; Qb = 1) at some value of R0 .
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 to accommodate the longest-wavelength unstable mode [6].In the incompressible system, where the injection rate is fixed at Qb = 1, Paterson [6] showed that this critical radius R0i satisfies Hence, R0i ∼ 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 Qb = 2 R0 Ṙ0 , where Ṙ0 ( t) is given by Eq. 7. The critical radius R0c in the compressible system must then satisfy 2Ca R0c where tc is the time of onset (i.e., R0 ( tc ) = R0,c ).To close Eq.( 21), we require R0 ( t), which we calculate numerically from the axisymmetric model.Hence, onset in the compressible system depends on both Ca and C, with the latter dictating the unsteady injection rate.
FIG. 10: Onset radius R0c in the compressible system [Eq.( 21)], normalised by the onset radius R0i in the incompressible system [Eq.(20)] at the same Ca, as functions of Ca for C ranging from 0.01 to 10 (increasing in the direction of the arrows).Red and blue curves, respectively, indicate C > 1 and C < 1.The dashed black line denotes R0c = R0i .Inset: Close-up of the low-Ca region where compressibility can weakly promote onset ( R0i < R0c ), according the predictions of linear stability analysis.
In Figure 10, we compare the linear stability predictions of onset in the compressible and incompressible systems by plotting R0c / R0i as a function of Ca for varying C. We determine R0c numerically, while R0i is given analytically by Eq. ( 20).We observe a substantial delay in the onset of the instability, indicated by R0c / R0i > 1, when Ca ≳ 10 2 .At the highest Ca and 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 Qb ≪ 1 [see Fig 5(a)] associated with compression-driven displacement.The effect is amplified with increasing C, which leads to lower initial Qb , and with increasing Ca, which corresponds to smaller onset radii in the incompressible system.Hence, if onset is predicted at small radii R0 ≪ 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 R0c / R0i at Ca ≲ 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 R0c / R0i < 1.This promoting effect at low Ca is due to the relatively high injection rates ( Qb > 1) at later times or larger R0 , 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 R0 ⪅ 1, and it is therefore unclear whether the effect on viscous fingering would be visible before breakout.

E. 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 where L and A are, respectively, the length of the interface and the area it encloses.For a circular interface, I = 1.Hence, any deviation from unity indicates some perturbation away from the axisymmetric base state assumed in the linear stability analysis [ § V D].FIG.12: (a) Onset radius R * at which I first exceeds 1.1 (c.f.Fig. 11) and (b) the near-breakout isoperimetric ratio I( t0.9 ) measured at t0.9 = t(max( R) = 0.9) as functions of C for a range of Ca.Symbols and colours correspond to Ca = 2.61 × 10 3 , 5.21 × 10 3 , 1.04 × 10 4 , and 2.08 × 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).
correspond to the top row of Figures 3 and 4, respectively.The isoperimetric ratio 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 I ≈ 1 occurs at larger radii [max( R)] for larger C (arrows) in both simulations and experiments, consistent with a delayed onset.The growth of I with max( 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 C impacts the onset of fingering, we define the radius at onset R * in our experiments and simulations as being max( R) at the last recorded instant (video frame or time step) for which I < 1.1.The measured and computed values of R * as functions of 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 R * tends to increase with increasing 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 D, confirming the prediction that increasing C delays the onset of fingering.Moreover, varying Ca or 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 FIG.13: (a) Volume of air in the flow cell at near-breakout, Vb ( t0.9 ), and (b) overall average injection rate up to near-breakout, Vb ( t0.9 )/ t0.9 , as functions of C for a range of Ca.Note that Vb ( t0.9 ) is also approximately equal to the near-breakout volume of liquid expelled since Vb (0) ≪ Vb ( t0.9 ).Symbols and colours correspond to Ca = 2.61 × 10 3 , 5.21 × 10 3 , 1.04 × 10 4 , and 2.08 × 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).
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 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 I( t0.9 ).That is, the isoperimetric ratio at the near-breakout time t0.9 = t(max( R) = 0.9) (see § III B).We plot I( t0.9 ) as a function of C for different values of Ca in Fig. 12(b).In simulations, we find that I( t0.9 ) tends to the incompressible case (dashed lines) as C → 0, and the severity of the fingering pattern at t0.9 decreases with increasing 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 B, 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 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 C. The near-breakout volume of air in the cell, Vb ( t0.9 ), which is approximately equal to the near-breakout volume of liquid expelled since Vb (0) ≪ Vb ( t0.9 ), is close to the incompressible value and essentially independent of C for C ≲ 0.1 [Fig.13(a)].Recall that the same is true of the breakout time tf [Fig.8(a)].Thus, the average injection rate up to near-breakout (i.e., the average value of Qb during the interval 0 ≤ t ≤ t0.9 ), which is given by Vb ( t0.9 )/ t0.9 , is approximately unity for C ≲ 0.1 [Fig.13(b)].That is, the average injection rate is approximately equal to the nominal injection rate for C ≲ 0.1, despite the non-trivial time evolution of the instantaneous injection rate Qb (t).Nonetheless, these smaller values of C lead to noticeably delayed onset and reduced intensity of fingering [Fig.12].Both Vb ( t0.9 ) and I( t0.9 ) increase substantially as C increases further (C ≳ 0.1) [Figs.12b and 13b], but Vb ( t0.9 )/ t0.9 decreases substantially [Fig.13].Thus, weak compressibility (C ≲ 0.1) delays onset and mitigates fingering by introducing a time-varying Qb ( t) while roughly preserving the average injection rate, whereas stronger compressibility (C ≳ 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 Vb ( t0.9 ) in experiments than in simulations [Fig.13a], which is consistent with the systematically larger values of I( t0.9 ) [Figs.12b] and further suggests that the experiments are subject to more severe fingering than the simulations at the same values of Ca and 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 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 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 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-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 C = 1, which dictates the transition between quasi-steady and burst-like dynamics.In our axisymmetric model, C plays a directly analogous role, with the trivial solutions of the system vanishing for R0 ≤ 1 at C = 1.However, due to the evolving base state (increasing radius) of the axisymmetric model, we identified an additional sensitivity to the initial radius 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 D, 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 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 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 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 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) [20,54].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 C ≲ 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 C ≳ 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 C ≈ 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 stateof-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 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 ε (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 [22,24,26,28,59] 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.On the one hand, when C ∝ QV res (0) is smaller, the difference between adiabatic and isothermal models is minimal [Fig.A1(a)] because the pressure ∆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 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 C is larger, the two models predict significantly different Vb for the same ∆p g [Fig.14(b)].Because Fo ≳ 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 Vb from the isothermal and adiabatic models against the normalised projected volume Â = Ab/(πR 2 c 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 Fo.

FIG. 4 :
FIG. 4: Example numerical solutions for Ca increasing top to bottom.Column one shows solutions of the incompressible model, while columns two to eight show numerical solutions of the compressible model Eqs.(3a)-(3e) with C increasing left to right.Profiles are plotted in equal time intervals of t 0.9 /10.Simulations with shaded background correspond to the same parameter values as the experiments shown in Fig. 3.

FIG. 7 :
FIG. 7: Phase-space representation of the axisymmetric model [Eq.(7)], in terms of injection rate Qb and interface radius R0 , for C = 0.25.The local time derivative of injection rate Qb is indicated by the colormap (see colorbar).Stable and unstable trivial solutions Q± are shown as solid and dashed black lines, respectively.The R0 and Qb coordinates of the saddle-node bifurcation point are indicated by dot-dashed and dotted lines, respectively.A solution of the axisymmetric model for C = 0.25 and R = 0.01 is plotted as a thick blue curve.Regions of the phasespace for which Qb is positive or negative are shaded orange or purple, respectively (see colorbar).When λ(∆p g ) > λ(ω), the interface accelerates, Qb > 0. This occurs in three distinct subregions (see inset): when Qb < 1 because the gas is compressing while resistance decreases (∆ ṗg > 0 and ω < 0); when 1 < Qb < Q− because the gas is decompressing slowly (∆ ṗg < 0 is small); and when Qb > Q+ because the gas is over-compressed (∆p g is large).When λ(∆p g ) < λ(ω), the interface decelerates, Qb < 0. This occurs only in the region Q− < Qb < Q+ , where the gas decompresses quickly (∆ ṗg < 0 and |∆ ṗg | is large).

FIG. 8 :
FIG. 8: (a) Breakout time at which the interface first reaches the rim of the flow cell and (b) corresponding breakout gauge gas pressure ∆p g ( tf ), both computed from numerical solution of Eq. (7) for R = 0.025 (black lines) and from the full numerical simulations (coloured curves and symbols) for Ca ∈ {2.61 × 10 3 , 5.21 × 10 3 , 1.04 × 10 4 , 2.08 × 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).The incompressible results are shown for reference (horizontal dashed lines).

FIG. 14 :Figure 14
FIG.14:The normalised total experimental gas volume Vg as a function of normalised time t for isothermal (solid blue; 13) and adiabatic (dotted red; A1) models of gas compression.Dashed black curves are the normalised projected volume Â( t) of the fingering pattern, which is an upper bound on Vg .Panels (a) and (b) show data from experiments performed at the smallest and largest compressibility numbers studied, C = 0.018 (Q = 1.25 mL/min, V g (0) = 25 mL) and C = 1.13 (Q = 10 mL/min, V g (0) = 200 mL), respectively.

FIG. 15 :
FIG. 15: Experimental reproducibility.(a, b) Interface evolution for experimental repetitions at Q = 2.5 mL/min and (left to right) V res (0) ∈ {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 (c) volume, (d) gauge gas pressure ∆p g , and (e) isoperimetric ratio as functions of non-dimensional time for the experiments shown in (a) and (b), plotted as solid red and dotted blue curves, respectively.Arrows indicate increasing V g (0).