Peeling fingers in an elastic Hele-Shaw channel

,


I. INTRODUCTION
Thin-film flows confined by elastic boundaries occur in a plethora of applications related to industrial, geophysical and biological processes [35], e.g., from roll coating [4] to spleenic filtration of red blood cells [36].In microfluidics, there is considerable interest in incorporating thin membranes and other soft boundaries into devices to harness flow-induced functionality and enable passive flow control [2,8,18,43].However, many practical microfluidic flows are multiphase, and are thus jointly influenced by elastic boundaries and moving fluid interfaces, which each introduce nonlinearities into otherwise linear flows in the Stokes limit.In this paper, we use a combined approach of experiments and numerical modelling to map out remarkably complex dynamics in a conceptually simple example of such a flow, namely a two-phase displacement flow in a Hele-Shaw channel where the top wall is an elastic sheet.
Two-phase displacement in a rigid Hele-Shaw channel, a channel whose width is much greater than its depth, is an exemplar of complex interfacial dynamics [5,9,29].When a less viscous fluid, usually air, invades a more viscous fluid, the Saffman-Taylor viscous fingering instability leads to the steady propagation of a single, centred finger [42], which can be captured accurately with numerical models [34], provided the liquid films left behind the advancing finger tip are taken into account [44].Although this finger is linearly stable [3] for all dimensionless finger speeds, quantified by a capillary number , finger instabilities are readily observed under finite-amplitude perturbations imposed either via the geometry [9] or via the fluid, e.g., by suspending particles [7].
In the rigid configuration, injected air displaces resident fluid along the length of the channel.In contrast, the injection of air into a liquid-filled compliant channel inflates its elastic top wall, imposing a reopening channel profile which localises fluid redistribution to a wedge-shaped region ahead of the advancing interface [32].
This results in the peeling of the elastic top sheet from the rigid bottom wall at an angle which is set by the fluid-structure interaction [37][38][39].Related peeling scenarios arise in the context of pulmonary airway reopening [20,27], where benchtop models have characterised the steady propagation of an air finger into two-dimensional elastic channels under axial tension [16,17,30] and collapsed elastic tubes [22,31].For moderate levels of tube collapse where the top elastic wall does not contact the bottom boundary, reopening takes place either via steady peeling modes, where an increase in pressure drives faster fingers, or steady pushing modes, where the converse is true, i.e., an increase in pressure drives slower fingers.Transition between these two modes occurs at a critical and a yield pressure difference (bubble pressure relative to the pressure in the collapsed channel), which must be exceeded in order for a finger to propagate.The less intuitive pushing mode ( < ), which has been found to be unstable in a two-dimensional model [21], follows from the fact that as is reduced, the elastic channel must expand indefinitely to accommodate redistribution of a finite volume of fluid within the ever-thinning films on its walls [22].
It is well established that in a radial Hele-Shaw cell of uniform depth, an elastic top wall leads to the suppression of viscous fingering instabilities [6,32].This is because the interface advances into a convergent channel, where both viscous and surface tension forces act to stabilise interface perturbations [1,37,40].The situation is more complex for finger propagation in our rigid Hele-Shaw channel where the top boundary is an elastic sheet.This is because the presence of side walls and initial channel collapse lead to a cross-section of non-uniform depth, with maximum constriction in the centre of the channel.The reopening mechanics of this system were first explored by [13].Remarkably, they found that unsteady finger propagation characterised by complex pattern formation at the finger tip, was supported for sufficiently large across a broad range of initial levels of collapse.For relatively large initial collapse, [10] identified two apparently persistent modes of finger propagation amongst a host of transient dynamics, which were mediated over a region of bistability by unstable pushing behaviour.These distinct reopening modes were shown to be dominated by viscous and elastic forces, respectively.The elastic reopening mode was associated with small-amplitude fingering perturbations which yielded intricate interfacial patterns referred to as 'feathered' modes.However, the nature of these feathered modes remains elusive.
A modelling framework for fluid-structure interaction flows was proposed by [14], based on a depth-averaged model to describe the propagation of an air finger into a collapsed elasto-rigid channel.They showed that the presence of the elastic wall can lead to interaction between solution branches that are isolated in the rigid channel, thus altering their stability and potentially leading to complex dynamics at higher levels of initial collapse.However, the model was unable to predict the myriad of exotic fingering patterns found experimentally [10,13].In this paper, we enable the simulation of intricate interfacial patterns by extending the depth-averaged model introduced by [14] to include an artificial disjoining pressure to prevent self-intersection of the interface.We systematically investigate the transition to feathered modes across a range of levels of collapse.Remarkably, we find that they arise after long oscillatory transients which resemble our experimental observations at lower capillary numbers.Furthermore, both experiment and model indicate that the small scale indentations which develop and advect around the finger tip are refined through tip-splitting as the finger propagates.This implies that the feathered modes of propagation are in fact continually evolving, with no evidence that these disordered states converge to steady or periodic states.[14] also demonstrated that a thin-film model is fundamental to capturing experimental results both qualitatively and quantitatively.In this paper, we refine their thin-film model to improve the prediction of finger speed as a function of injection flow rate, but we also find that the dynamics are not sensitive to the exact choice of thin-film model.
The paper is organised as follows.We recall the experimental methodology in §II and highlight the essential new features of the numerical model in §III.Results are presented in §IV.We introduce a measure of the dimensionless finger speed based on viscous dissipation within the fluid ahead of the finger tip in §IV A, which enables us to identify a threshold value of beyond which viscous forces are superseded by elastic effects.Quantitative prediction of this transition between 'viscous' and 'elastic' reopening regimes across levels of collapse establishes the fidelity of the numerical model.We show in §IV B that in the viscous regime, we recover the non-monotonic dependence on of the finger pressure, which is characteristic of benchtop models of airway reopening.In §IV C, we explore the elastic regime and the transition to feathered states as a function of .We capture feathered states numerically for the first time in §IV C 1. A detailed analysis of the steady bifurcation structure in §IV C 2 reveals that steady numerical solutions match the bounding envelopes of unsteady fingers in the elastic regime and, thus, satisfactorily predict the bubble pressure.However, the steady bifurcation diagram does not inform the transition to unsteady dynamics, which appears strongly nonlinear.

II. EXPERIMENTAL SET-UP
We performed experiments in a rigid Hele-Shaw channel topped with an elastic sheet, shown schematically in figure 1(a).This experimental set-up was previously described by [10].A channel of length 60 cm, width * = 30 ± 0.02 mm and depth * 0 = 1.05 ± 0.01 mm was precision-milled into a block of Perspex, achieving a 10 m roughness along the base.For the top boundary, we used a Latex sheet (Supatex) with Young's modulus * = 1.44 ± 0.05 MPa, Poisson's ratio = 0.5 and thickness ℎ * = 0.46 ± 0.01 mm.A uniform pre-stress was imposed on the elastic sheet, directed across the width of the channel, by hanging evenly distributed weights (totalling 3.03 kg) along one long edge of the sheet.This pre-stress was maintained by clamping the elastic sheet in place using an aluminium frame secured with 11 evenly-spaced G-clamps along each long edge and a single bolt along each short edge.The constitutive relation between transmural pressure difference across the elastic sheet and level of collapse, which was measured in the experimental channel under static conditions, matches numerical simulations performed using a pre-stress that is 30% of the value predicted from the sheet dimensions and applied weight (see appendix A).This loss of applied pre-stress is expected, due to unavoidable slippage of the elastic sheet during the clamping procedure.However, the reproducibility and consistency of the experimental results suggest that the pre-stress was uniform and remained constant once the elastic sheet was clamped in place.
To set the initial conditions of each experiment, the channel was filled with silicone oil of dynamic viscosity * = 0.099 Pa.s, surface tension * = 21 mN/m and density * = 973 kg/m 3 at laboratory temperature * = 21 ± 1 • C. The inlet was then closed and oil was allowed to drain from the outlet of the channel, via a length of flexible tubing open to the atmosphere at one end.By adjusting the height of the outlet of the tubing relative to the channel we could set the hydrostatic pressure difference between the channel and the atmosphere.At equilibrium, this hydrostatic pressure difference is equal to * trans , the transmural pressure difference across the elastic sheet, which acts to deform the sheet.Hence, we were able to control the extent to which the elastic sheet was initially collapsed.We measured the initial collapse using a laser sheet projected to a line across the width of the channel and imaged this line at an oblique angle using a camera of resolution 22.9 ± 0.02 pixels/mm (the same camera was used for visualising sheet deformation during the reopening experiment, see below).The resulting profiles of the sheet in terms of the local depth * of the channel as a function of the lateral coordinate * 2 are shown in figure 1(b).We quantify the extent of initial collapse with the parameter , which is defined as the ratio of cross-sectional area under the sheet to the uncollapsed area * * 0 .Hence, for an uncollapsed channel = 1.By = 0.36, the elastic sheet sags sufficiently to contact the centreline of the bottom wall, which we refer to as opposite wall contact.We studied channels collapsed within the range 0.43 ± 0.02 ≤ ≤ 0.95 ± 0.02, as detailed in figure 1(b).
The channel was reopened by injecting air at a constant volumetric flow rate * , which varied within the range 5 ≤ * ≤ 330 ml/min, resulting in a long continuous finger of air.The finger propagated along the length of the channel, displacing oil and parting the channel walls.Flow was imposed using a syringe pump (KD Scientific) fitted with Gastight syringes (Hamilton), with a small precursor bubble (< 1 mL) injected immediately before each experiment to set the initial conditions.The pressure * b of the air finger was measured with a differential pressure sensor (Honeywell, ± 5" H2O) with ±1 Pa resolution.Over a 250 mm length of the channel, which was chosen as the region of interest (ROI) for our measurements, we recorded constant air-bubble pressure traces to within a typical tolerance of 10 Pa.
The air finger was imaged from above by a second camera (4.8 ± 0.1 pixels mm −1 ) recording images at fixed rates of 0.33-30 frames per second (f.p.s.) depending on * and .The channel was back-lit by LED lights diffused through opalescent acrylic.An example of a top-view image is shown in the upper panel of figure 1(c).The tip of the air finger was located in each frame using image analysis routines (MATLAB 2016a), which allowed us to calculate the instantaneous axial speed * f of the finger.To illustrate the time-evolution of the finger over an experiment, we will refer to composite images of the kind shown in the lower panel of figure 1(c), which are generated by overlaying background-subtracted images from successive frames of an experiment.In figure 1(c) the finger shape is steady over the entire ROI and, since the interfaces are equally spaced at fixed time intervals, the finger advances at a constant speed.
The deflection of the elastic sheet during reopening was measured by recording images of the laser sheet (at 40-160 f.p.s.), which was located near the end of the ROI, see figure 1(c).Measurements of the height * of the elastic sheet mid-way between the walls of the channel along with knowledge of the finger speed * f allowed us to reconstruct the profile of the elastic sheet around the finger tip, as shown, e.g., in figure 1(d).Ahead of the finger tip (located at axial coordinate * 1 = 0) the channel is still collapsed, while behind it the channel is inflated.Hence, the air finger advances into a tapered region behind a wedge-shaped volume of oil.

III. MODEL
We extend the modelling framework developed by [14] which couples depth-averaged lubrication equations for the fluid flow to the Föppl-von Kármán plate equations with in-plane pre-stress for the elastic sheet.The original model was developed in a frame of reference moving with the axial velocity of the finger tip so that steady solutions represented fingers propagating steadily in the frame of the laboratory.This model was implemented numerically in oomph-lib [26] to calculate steady solutions, evaluate their linear stability and perform time simulations.Compared to the model of [14], the new model includes an additional artificial disjoining pressure, which has been implemented to prevent the unphysical self-intersection of air-finger interfaces when deep indentations develop during time-simulations, see §III B 2. Furthermore, the effect of different approximations for the liquid films separating the finger from the walls of the channel in comparison with experiments has been analyzed.

A. Governing equations
We define a frame of reference in Cartesian coordinates ( * We use the channel width, * , to non-dimensionalise the in-plane coordinates ( 1 , 2 ), and the undeformed channel height, * 0 , for the out-of-plane coordinate 3 .These define the aspect ratio = * / * 0 of the channel.The displacement in the elastic sheet is non-dimensionalised using the in-plane length scale and the fluid pressure is non-dimensionalised using P * = 12 * 2 /T * , where T * is the time scale defined as T * = * 2 * 0 / * .Finally, in-plane velocities are non-dimensionalised by U * = * /T * .Assuming incompressibility, we apply the Reynolds lubrication approximation to the Stokes equation, defined in the frame moving with instantaneous velocity f ( ) = * f ( )T * / * , which yields the depth-averaged governing equation for the fluid pressure : where we use the summation convention and the index , and all subsequent Greek indices, takes the values 1 and 2. The unknown speed f ( ) is determined by enforcing that the maximum 1 coordinate of the interface (i.e., the 1 coordinate of the finger tip 1,tip ) is set to zero.We use the Föppl-von Kármán equations [33] 2 2 as the governing equations to determine the in-plane and out-of-plane displacements of the elastic sheet, ( 1 , 2 ) and , respectively.The pressure load * on the elastic sheet is non-dimensionalised using the bending modulus * = * ℎ * 3 12(1− 2 ) so that = * * 3 / * .The parameter = 12(1 − 2 ) * ℎ * 2 describes the relative importance of the in-plane and bending stresses.Finally, the in-plane components of the stress tensor, are 11 = where the in-plane strain is given by and a non-zero pre-stress component (0)  22 is introduced to mimic the clamping procedure performed in the experiment which tensions the sheet in the * 2 direction.The pre-stress components 11 and (0) 12 are fixed as zero.However, as mentioned in §II the experimental value of (0) 22 is not known accurately, so, following [14], we estimate it by matching the constitutive relation of the experimental channel to numerical simulations.Details of this procedure are presented in appendix A.
The equations governing the fluid pressure (1) and elastic sheet deformation (2) are coupled in two distinct ways.First, the height of the channel is determined by the out-of-plane displacement of the elastic sheet: Secondly, the pressure load on the elastic sheet depends on the pressure of the fluid and the pressure b of the air finger: where the non-dimensional fluid-structure interaction parameter provides a measure of the typical viscous stresses in the fluid relative to the bending stress of the elastic sheet.As I → 0 the elastic sheet becomes rigid and the governing equations ( 1) and ( 2) decouple.

B. Boundary conditions
At the side walls of the channel, the boundary conditions are non-penetration of fluid and a clamped elastic sheet: The numerical domain is truncated in the axial direction at the upstream ( 1 = − up ) and downstream ( 1 = down ) ends, see figure 2(a).We set up = 10 and down = 15, and checked that the computed solutions are not affected by increasing the length of the domain.Following [22], at these truncated boundaries, we impose the conditions These mean that far away from the finger tip, all disturbances should decay.We allow the upstream pressure gradient to be an unknown, which we determine by imposing that the fluid flux at the downstream boundary is equal to the flux at 1 → ∞: where is the initial level of collapse defined in §II.

Effect of fluid films left behind the advancing interface
Following [38] and [14], we incorporate into the boundary conditions at the air-finger interface the effects of fluid films left behind the advancing interface.The kinematic boundary condition is given by where R( , ) is the position of the advancing air-fluid interface in the moving frame, parameterized by the coordinate , and n is the in-plane outer unit normal vector to the interface, see figure 2. The dynamic boundary condition is given by where is the in-plane curvature of the air-finger interface and = * * f / * is the capillary number.Note that the capillary number is based on the instantaneous velocity of the finger tip.This means that, for unsteadily propagating fingers is timedependent.The functions 1 ( ) and 2 ( ) incorporate the effect of the fluid films into the model.We use the expressions and which [38] introduced and validated for use in elastic cells.The effects of fluid films are removed if we set 1 ( ) = 0 and 2 ( ) = 1.[38] assumed that the thickness of the fluid film is set at the finger tip position Π tip = ( 1,tip , 2,tip ), where it is formed.This means that, in their model, the fluid film has a uniform thickness 1 ( ) tip at every point ( 1 , 2 ) of the air finger, where tip = ( 1,tip , 2,tip ).This choice differs from the assumption made by [14], where the thickness of the fluid film was a constant proportion of the channel height 1 ( ) ( 1 , 2 ) at each point ( 1 , 2 ) of the air finger.We refer to these assumptions as uniform and non-uniform film thickness models.In this paper, we follow [14] and assume a non-uniform film thickness model.In appendix B we compare both assumptions with experiments and, in addition, show that the relationship between bubble pressure and capillary number predicted by the model remains the same under both assumptions.The imposed flow rate necessary to reach a given capillary number, however, depends on the assumption made for the film thickness.In fact, both assumptions are simplifications of the three-dimensional fluid rearrangement that takes place within the films; and three-dimensional Stokes simulations would be needed to accurately capture the distribution of fluid films required for detailed quantitative agreement with experiment.Moreover, we find that the behaviour of the fluid in these films does not affect the overall qualitative dynamics.
2. Disjoining pressure [14] found that time simulations at high and low produced unsteady fingers with small scale indentations originating near the finger tip.These indentations could grow into clefts sufficiently deep and narrow that they led to the self-intersection of the air-finger interface, thus terminating the simulation.However, self-intersection of the air-finger interface was never observed in the experiments presented by [13], [10] and in this paper.It is most likely that the self-intersection in the model is a result of the approximations made when simplifying the three-dimensional liquid-film dynamics.In order to allow time simulations to continue beyond the point of self-intersection, rather than moving to three-dimensional calculations, we prevent self-intersection by adding an artificial repulsive disjoining pressure d to our model.We choose the form based on the disjoining pressure in thin films introduced by [11].The constant gives the strength of the repulsive interaction and is the length scale of the interaction.The interface-interface distance , see depiction in figure 3(a), was calculated as the distance between a pair of points on the interface, one on each side of the indentation.Each interface point was paired with the point with the shortest distance from it and with the additional constraint that the unit normal vectors to the interface at these two points have a negative inner product.This constraint ensures that the points are always localized at opposite sides of the indentation.Once a pair is assigned for each point at the interface, we calculate the local value of the disjoining pressure d for that point.In order to identify the pairs of points located at indentations, we map all pairs in the entire finger interface.It is possible to show that for the pairs of points that are not on an indentation (e.g., a pair of points located far behind the finger tip, on opposite sides of the finger width), the distance is always of the order of the finger width (between 0.4 and 0.7), while for pairs that are on an indentation, the distance is 0.05 or smaller.For sufficiently large values of , d becomes negligible and, in order to speed up the calculations by reducing the length of the interface that must be scanned, we set the value of d to zero for points separated by a distance greater than a cutoff value of 0.1.For our chosen values of and , see below, < 10 −6 when > 0.1.The overall effect of this implementation of disjoining pressure is a short ranged repulsive interaction that only acts to prevent self-intersection on interface indentations.
The disjoining pressure was added to the dynamic boundary condition (12), resulting in a new condition We used trial and error to select the parameter values = 10 −3 and = 10 −2 , so that the disjoining pressure prevents self-intersection of the interface during time-dependent calculations but only affects the dynamics when the interface is on the verge of self-intersection.For steady state simulations, we set d = 0.
Figure 3 shows a comparison between time-simulations without (figure 3(b)) and with (figure 3(c)) the addition of a disjoining pressure.Both simulations are shown for the same time step and simulated with the same parameters and initial conditions.In the time-step following the snapshot of figure 3(b), the interface self-intersects, thus terminating the numerical simulation, whereas this does not occur in figure 3(c).There, the addition of d reduces the depth of the indentation and increases the distance between the sides of the indentation.This brings unsteady finger dynamics in the simulations closer to the experiments, where indentations reduce in depth as they are advected to the side of the finger before eventually disappearing; see §IV C 1.However, 3D Stokes simulations would be necessary to fully capture the interface-interface interaction in the region of indentations.In the case of deep indentations, the separation * can be as small as the unperturbed channel gap * 0 , making the lubrication approximation invalid.

C. Numerical implementation
The numerical implementation of our model followed that of [14], with the addition of the disjoining pressure to the fluid dynamic boundary condition and the possibility to choose between fluid film correction of uniform and non-uniform thickness.
In order to calculate steadily propagating modes (travelling-wave solutions), we set all time-derivatives in the governing equations and boundary conditions to zero.Additionally, we imposed fixed values of the initial level of collapse and another global variable which could be either the bubble pressure b , capillary number or flow rate * .Different choices were made for the second controlled variable according to the parameter continuation between solutions branches in the bifurcation diagram to be calculated.Once we had steady solutions, we analysed the linear stability of these solutions at fixed * and for consistency with the experiment.The solution of the discrete generalised eigenproblem was obtained via the Anasazi solver from Trilinos [28].We computed the 60 most unstable eigenvalues.
Lastly, in order to assess the non-linear stability of the steady solutions, we conducted time simulations for fixed values of the initial level of collapse and flow rate.We used a steady solution as initial condition to which we applied a time-decaying perturbation to the pressure jump across the interface: where the perturbation takes the form This perturbation creates a dimple on a length scale centred on an interfacial point with coordinate 2 = 0 .This artificial extra pressure is zero at = 0, reaches the maximum 0 and decays to zero for ≫ .We have chosen this particular form of perturbation in order to be able to apply it consistently for all values of parameters used.We set a non-zero value for 0 , so we prescribe a perturbation asymmetric about the centreline of the channel, thus avoiding restricting time-simulations to symmetric fingers.In our simulations, we used the following parameters: = 0.035, 0 = 0.05, = 0.04 and 0 = 0.2 b .

IV. RESULTS
We report experiments and numerical simulations of bubble propagation in a collapsed channel, see figure 1, for different values of and the imposed air flow rate * .Opposite wall contact is avoided by keeping the level of collapse above the value 0.36.We present results in terms of the capillary number, 0 < < 1.3, which provides a measure of the dimensionless finger velocity and is not sensitive to the choice of film model as discussed in §III B 1. The choice of experimental parameters listed in §II means that the fluid-structure interaction parameter I varies only with * in the range 0 < I < 1.5 × 10 4 , and that the other non-dimensional parameters are fixed at = 28.6 and ∼ 40, 000.

A. Comparison with the flow into a rigid wedge
We begin by examining the reopening region ahead of the finger.We define a flow metric which we use to establish the overall fidelity of the model by comparing experimental and numerical reopening flows across different and levels of initial collapse.When air is injected into a fluid-filled, collapsed elastic channel, the cross-sectional area of the collapsed channel expands and the fluid within it redistributes to accommodate finger propagation.This process takes place in a reopening region ahead of the finger tip where the local fluid-structure interaction flow sets the mode of finger propagation [10,32].If we assume that the reopening is dominated by viscous stresses, then we can approximate the process as flow into a rigid wedge (RW).This geometry is illustrated in figure 4(a), which shows a schematic diagram of the vertical mid-plane of the channel ( * 2 = 0) where the wedge-like reopening region has angle , length Δ * and pressure drop Δ * .At low , flow into a compliant channel is closely captured by the flow into a rigid convergent laterally unbounded channel of pressure gradient Δ * /Δ * , where the wedge angle is set by fluid-structure interaction [16,30,37].This approximation is completely 2D, and since there is no in-plane curvature, only the transverse curvature contributes to the pressure jump across the interface.In this rigid wedge, the capillary number corresponding to the dimensionless tip speed is directly proportional to the pressure gradient across the wedge and, following the depth-averaged approach, can thus be defined as Here, * ( 1,tip , 0) = * 0 ( 1,tip , 0) is the dimensional height of the elastic sheet at the centreline of the channel and at the 1 coordinate of the finger tip, the pressure drop is approximated by where * b is the air pressure in the finger and * trans the pressure in the collapsed channel far ahead of the finger, and the length of the wedge Δ * can be estimated as where * (∞, 0) = * 0 (∞, 0) is the dimensional height of the membrane at the centreline of the channel, beyond the wedge region (see figure 4(a)).The geometric quantities * ( 1,tip , 0), * (∞, 0) and were measured while the finger crossed the ROI, and remained constant even for unsteady finger propagation.
Figure 4(b) shows a plot of RW , which is estimated in our elastic channel based on the pressure gradient across the wedge, as a function of the actual capillary number based on the finger speed.In the experiment, the finger speed * f could exhibit measurable fluctuations depending on parameters and thus, we use its time-averaged value * f , over the time interval during which the finger-tip propagates in the ROI, to calculate a time-averaged capillary number = * * f / * .For the smallest values of investigated ( < 0.1), we find RW ≃ as expected in a laterally unbounded rigid wedge.For low levels of collapse ≥ 0.8, RW is approximately proportional to over the range of investigated.The reduced slope of the RW curve, compared to that of rigid wedge, is due to changes in the finger geometry in the tip region.The fact that there is linear relationship indicates that the fluid redistribution associated with reopening is primarily shaped by viscous and capillary forces, so we refer to this regime as viscous reopening.This viscous reopening regime also occurs for larger initial collapse, with the datasets for = 0.60, 0.53 and 0.43 exhibiting quasi-linear behaviour as a function of for small values of .As increases, the gradient of the RW curve for = 0.60 progressively decreases.This trend is enhanced for = 0.53, in which case RW becomes approximately constant at a threshold value, ≈ 0.17, indicating saturation of the pressure gradient over the fluid wedge, and thus of the viscous stresses, while the velocity of the finger continues to increase.The pressure gradient saturates because both the elastic wall and air finger have reached their limiting geometric configurations.Therefore, the fluid wedge no longer changes its shape and the pressure drop across it remains constant on the viscous scale, as also found in three-dimensional elastic tubes by [22] for sufficiently high propagation speeds.In fact, for = 0.43, RW drops to values close to zero for > 0.1 (red points in figure 4(b)), which suggests that the influence of viscous stresses becomes negligible for high levels of collapse approaching the point of opposite wall contact.This is because the elastic sheet steepens significantly in the reopening region, leaving it occupied mostly by the peeling finger.In this configuration, the value of * ( 1,tip , 0) becomes so small that the driving pressure difference ( * trans − * b ) is counteracted by a large capillary pressure 2 * * ( 1,tip ,0) .The result is a small estimated pressure gradient Δ * Δ * and therefore, a small RW .Thus, the finger is primarily shaped by elastic and capillary forces [10,12], so we refer to this regime as elastic reopening.
The comparison between experiments and steady numerical simulations is detailed with individual plots of each level of collapse in figure 4(c).Since these are steady numerical simulations, is used as .A comparison is not shown for = 0.43, because of difficulties in resolving the simulations for levels of collapse very close to opposite wall contact, when the fluid layer in the channel becomes very thin.However, the numerical model (solid lines) shows remarkable quantitative agreement with the experimental data (circles) for levels of collapse 0.53 ≤ ≤ 0.95.This indicates that our lubrication-based fluid-structure interaction model captures the key physics underlying finger propagation in the experiment.We find an initialcollapse-dependent threshold value of that divides viscous reopening from elastic reopening.Viscous reopening occurs at low and is the only observed behaviour at low levels of initial collapse.Elastic reopening occurs at high when the initial collapse is sufficiently large.We now proceed to discuss the modes of finger propagation associated with these different reopening regimes.

For
= 0.95 and = 0.82, all experiments gave steadily propagating fingers, each of approximately constant capillary number, consistent with previous studies of airway reopening [13,20,24].Figure 5 shows sequences of time-lapse images from experiments for = 0.95.Steady propagation is indicated by the fact that the shape of the finger does not vary as it propagates and that the consecutive interfaces are uniformly spaced.As the flow rate is increased from * = 5 ml/min (panel 1) to 50 ml/min (panel 2), the symmetric finger narrows and its tip curvature increases.This trend is reversed upon further increase of the flow rate for * ≥ 150 ml/min (panels 3-5) with the finger widening and reducing its tip curvature.The red interfaces in each panel are finger profiles from steady numerical calculations with the same capillary number as the experiments.They accurately capture the finger profile in panels (1-3) but differences appear in panels (4)(5), where the numerical solution is asymmetric about the centreline of the channel, in contrast to the experimental fingers which remain symmetric.
The non-monotonic variation of the finger width as a function of is plotted in figure 6(a), with good agreement between experiments (red symbols) and steady simulations (black lines).The finger width only decreases with increasing for the smallest values investigated, < 0.1.This is the limit where tends to in §IV A and a decreasing finger width concurs with displacement flows in rigid Hele-Shaw channels [44].The steady simulations undergo a pitchfork bifurcation near = 0.55, where symmetry about the centreline of the channel is lost.This bifurcation is responsible for the morphological change seen between panels (3) and (4) of figure 5.In figure 6(b), the steady simulations (solid line) capture the overall increase of the experimental bubble pressure (red circles) with increasing .The turning point at ≃ 0.15 in the numerical curve is a signature of compliant channel reopening, which indicates a transition from pushing to peeling fingers [10,16,22].In the viscous pushing regime, the * b curve has a negative slope and a large volume of fluid is displaced by the finger tip.In the viscous peeling regime, the * b slope becomes positive and very little fluid is pushed by the finger tip.In both these regimes, however, RW is proportional to , indicating the dominant balance is between viscous and capillary stresses.In contrast to the numerical results, the experimental data does not indicate a turning point at low , but the bubble pressure fluctuates significantly for the smallest value of = 0.01, with variations of ±20% over the ROI, despite an approximately constant velocity.This may be attributed to the presence of gravity forces in the experiment [23].The increasing bubble pressure of peeling fingers with increasing means that the reopening height of the channel increases and thus, volume conservation requires the fluid behind the finger tip to occupy a smaller fraction of the channel width [13].Hence, the finger width increases with as shown in figure 6(a).We note that the variability of the capillary number over the ROI increases with (up to ±7%).This is consistent with thinner liquid films around the peeling finger, which makes the finger more sensitive to channel imperfections that can affect its velocity [10].The experimental data point with the highest in figure 6(b), is not in good agreement with the numerical prediction.We attribute this discrepancy to the fact that the experimental bubble pressure (500 ) is far outside the range of pressures (−200 < < +200 ) for which the constitutive channel law used in our model was calibrated, see appendix A for details of the calibration.The elastic sheets used in the experiments stiffen under inflation, meaning that in order to inflate the channel by the same amount higher pressures are required in the experiments than in the model; the level of channel inflation is set by conservation of mass.The stiffening increases nonlinearly, explaining why there is such a dramatic increase in the difference between experimental and numerical pressures at the highest value of .At = 0.82, we observe similar trends in both experiments and simulations; see Appendix C.However, the turning point in the numerical simulations is displaced to a 40% lower value of the capillary number.Moreover, a pitchfork bifurcation that leads to asymmetric fingers [14], indicated in figure 6 for = 0.95 by the emergence of a solution branch at ≃ 0.55, is displaced to smaller for = 0.82 and a symmetry-breaking bifurcation is also observed in the experiments.

C. Elastic reopening in strongly collapsed channels
We now turn to the elastic reopening regime which was identified in §IV A for increasing levels of collapse.Figure 7 shows experimental sequences of time-lapse images of finger propagation at = 0.36, for five values of the initial level of collapse 0.43 ≤ ≤ 0.78.Steadily propagating fingers are only observed in the least collapsed channel with = 0.78 (panel 1), which reopens viscously based on the findings of §IV A. By contrast, finger propagation is unsteady for the larger levels of collapse (panels 2-5) where a dominantly elastic reopening regime is expected.The unsteadiness is associated with the appearance of interface dimples on the finger tip, which can grow into clefts and advect around the curved tip as the finger propagates.They arise either intermittently or approximately periodically.For the largest initial collapse ( = 0.43), the interface indentations form small-scale corrugations on the finger tip which cause the characteristic feathered pattern previously highlighted by [10].
The disjoining pressure condition introduced in §III enables us to extend the numerical simulations benchmarked by [14]  to capture unsteady finger propagation prevalent in the elastic reopening regime.We focus on unsteady finger propagation at = 0.53 and compare numerical findings with experiments in §IV C 1 .We then compute the underlying bifurcation structure in §IV C 2 to gain insight into the transition from viscous to elastic reopening and unsteady propagation.

Multiple modes of finger propagation
Figure 8(a) shows the wide variety of finger propagation observed experimentally for = 0.53 and flow rates between 5 ml/min in panel 1 and 330 ml/min in panel 10 (0.01 ≤ ≤ 0.65).We classify the morphological structures seen in the experiments into four categories: (i) round-tipped, a finger propagating steadily with a symmetric curved front, depicted in panel 1 of figure 8(a); (ii) flat-tipped, a finger propagating steadily with an asymmetric flat front, depicted in panel 2 of figure 8(a); (iii) pointed, a finger propagating steadily with a symmetric triangular shaped front, depicted in panel 5 of figure 8(a); and (iv) feathered, a finger propagating unsteadily with small-scale perturbations continuously developing at the tip, depicted in panel 10 of figure 8(a).
For the lowest values of (panels 1 and 2), channel reopening is associated with the steady propagation of round-tipped symmetric and flat-tipped asymmetric fingers similar to those discussed in weakly collapsed channels in §IV B. For = 0.10 (panel 3), the asymmetric finger is prone to intermittent growth of tip perturbations, which is reminiscent of tip-splitting events in viscous fingering in rigid channels [9].The flat-tipped asymmetric state transitions to a pointed finger in panel 4 ( = 0.16).The appearance of the pointed finger, which is observed over the entire ROI in panel 5, coincides approximately with the saturation of in figure 4, which marks the transition to the elastic reopening regime.These pointed fingers were not observed within the range of flow rates investigated for = 0.60 (appendix C).Panels 6-10 in figure 8(a) all show unsteady finger propagation ( ≥ 0.24).In panel 6, dimples appear almost periodically on one side of the pointed finger tip and thus the finger propagation is asymmetric about the centreline of the channel.These periodic indentations grow into clefts in panel 7.By panel 8, the interface dimples form on both sides of the finger tip in alternation and the pattern of finger propagation appears mildly disordered.As discussed in §IV B, the finger width increases with increasing thus reducing the curvature of the finger tip.This leads to feathered modes of propagation in panels 9 and 10 because the weakly curved front of the finger tip promotes the appearance of small stubby fingers which form a regular corrugation [12].Variation of the characteristic width of these interfacial protrusions along the ROI could be due to imperfections in the channel but comparison between panels 9 and 10 suggests that their width decreases overall with increasing .The finger propagation is unsteady because the weakly-curved finger tip advects the small finger-like protrusions sideways from the centreline, thus broadening the central protrusions which in turn undergo tip-splitting.This results in a delicate cycle of small-scale stubby finger generation on the weakly curved front, where the length scale of the fingering decreases continuously as the finger propagates (panels 9, 10).Remarkably, the interface in these experiments is not prone to self-intersection despite its fine indentation, as discussed in §III B 2. This was also the case with the previously mentioned clefts in panels 7 and 8, which always remained sufficiently short and broad.
Figure 8(b) shows time simulations performed for the same values of mean capillary numbers as in the experiments in column (a), except for panels 9 and 10 where is smaller then in the corresponding panels in column (a).We made the choice of .53 and constant flow rate that increases from 5 ml/min in experiment 1 up to 330 ml/min in experiment 10.The time interval between the interfaces are, from 1 to 10: 1.17, 0.67, 0.40, 0.17, 0.17, 0.17, 0.17, 0.10, 0.17, 0.13 and 0.10 s, respectively.The third to last interface in experiment 6 is missing due to a camera fault.See also Supplementary Movies.(b) Time evolution of our time-dependent simulations at a fixed = 0.53 and constant flow rate that increases from 3 ml/min in simulation 1 up to 160.3 ml/min in simulation 10.The time interval between the interfaces is 0.08 in the non-dimensional scale, which results in 1.486, 0.205, 0.134, 0.078, 0.062, 0.052, 0.043, 0.038, 0.037, 0.028 s from 1 to 10, respectively.The initial steady solution in blue and the perturbed profile at in red are the first two finger profiles in each time-sequence.See also Supplementary Movies.taking smaller steps in for the simulations, so we could capture the change in the transient behaviours that precedes the onset of the feathered mode.Time-simulations were initialised with steady solutions (blue contours) which were similar to the finger envelopes observed experimentally: round-tipped symmetric fingers (panels 1, 3), flat-tipped asymmetric fingers (panel 4) and pointed fingers (panels 5-10); see §IV C 2. The distorted finger shapes after perturbation of the steady solution using equation ( 18) are shown with red contours at = = 0.04.In panels (1-7), steady modes of finger propagation are recovered over a distance of the order of a channel width.The transitions from symmetric round-tipped (panels 1, 2) to asymmetric flat-tipped (panels 3, 4) fingers, and to symmetric pointed fingers (panels 5-7), occur for similar values of as in the experiment.The absence of finger indentations in the simulations in panels (1-5) means that the disjoining pressure did not affect the evolution of these fingers.In the simulations in panels (6-7), the initial perturbation evolved into an indentation but interface self-intersection at this indentation was prevented by the disjoining pressure, thus allowing the simulations to proceed.
The first instance of unsteady finger propagation occurs at = 0.31 (panel 8) in the simulations and at = 0.24 (panel 6) in the experiments.For larger respective values of , both experiments and time-simulations only support unsteady finger propagation.For sufficiently large values of , regardless of initial transient, the finger develops into the feathered mode.Due to long transients, simulations in panels 8 and 9 only show the early and later stages of the finger propagation.Figures 8(c,d) show the full length of the numerical domain for these two panels.They indicate that the numerical model captures the unsteady dynamics observed in the experiment, where the pointed finger develops periodic perturbations (Figure 8(a), panels 6, 7) .In figure 8(c), asymmetric oscillations destabilize the symmetric pointed finger after a long transient.The perturbation develops a cleft close to the fingertip, which is advected and undergoes tip-splitting at approximately the same time.This process is repeated resulting in the feathered pattern.In figure 8(d), we observe a similar transition from pointed finger to feathered mode, but the transient oscillatory perturbation is symmetric instead.This transition from asymmetric to symmetric perturbations is seen in the experiments (Figure 8(a), panels 7, 8) with similar values of .However, the finite length of the experimental channel means that we cannot say for certain if the patterns seen in panels 7 and 8 are long term transients or periodic states.
Finally, at = 0.40 (figure 8(b), panel 10) the simulation depicts a finger that evolves from the symmetric pointed to the feathered state after a very brief transient.The time-dependent simulations of the feathered modes are interrupted when the width of the small-scale finger becomes comparable to the channel height * 0 .At this point the depth-averaged model is no longer reliable.As in the experiments, the typical lengthscale of the small-scale fingering seen in the simulations decreases as increases.However, the simulations feature small-scale indentations that are consistently deeper than in the experiments.Given that the typical indentation width is comparable to the unperturbed channel height * 0 , there will be three-dimensional effects in the experiments that cannot be captured by our model.At this stage, regardless of the calibration or the choice of functional form for the disjoining pressure, three-dimensional Stokes simulations are necessary to correctly capture the interfacial interaction between the small-scale fingers.
The propagation of the experimental feathered modes gives no indication that a periodic state is being reached, instead we can see that the length scale of the fingering pattern is always refining for the entire length of the channel.Similarly, the simulations of the feathered mode stop once the length scale of the fingering pattern reaches the length scale of the channel height.In figure 9 we gauge the apparent periodicity of the feathered modes by tracking the time evolution of the 2 component of a point Π = ( 1,tip − , 2 ( )) on the finger interface at a fixed axial distance = 0.5 behind the finger tip Π tip = ( 1,tip , 2,tip ).For both = 0.24 and = 0.28, after a transient due to the advected perturbations, 2 component of Π eventually reaches a steady value.For the feathered modes, figures 9(c)-(d), although the oscillations persists as long as the simulations last, they never reach a periodic state.This gives us an indication that the simulations point to the same result seen in the experiments, that there is no periodic state being reached.
As previously reported by [14], the model system exhibits a complex solution structure, with multiple co-existing steadily propagating states and numerous bifurcations.The turning point in the pressure that indicates a transition to pushing in the numerical simulations is at the lower end of the capillary number range investigated ( ≃ 0.03).The symbols of different colours distinguish the different finger types observed experimentally.As shown in figure 8, the steadily propagating round-tipped symmetric finger (solid green circle), flat-tipped asymmetric finger (solid red square) and pointed symmetric finger (solid blue triangle) occur in succession as increases.The steady numerical simulations capture their bubble pressures and ranges quantitatively.The insets show examples of experimental fingers overlaid with a steady numerical solution at the same .In insets 1, 2, and 3 where the finger propagation is steady, there is excellent agreement between simulations and experiments for finger width and morphology.The unsteady experimental fingers discussed in §IV C 1, which occur in the vicinity of transitions between modes and for > 0.27, are shown with open symbols.In the numerical simulations the solution branches for flat-tipped fingers (red), which include both symmetric and asymmetric solutions, and pointed symmetric fingers (blue) continue to values of the capillary number where only feathered fingers are observed in the experiment ( > 0.27).For up to 0.5, the bubble pressure of the experimental feathered modes lies between the red and blue curves but closer to the latter.As increases and complex eigenvalues with positive real parts , corresponding to instabilities, are indicated in brackets ( , ) for each solution branch.The relevant bifurcations are marked with red dots. 1 and 2 are pitchfork bifurcations, 1 -5 are Hopf bifurcations, 1 is a limit point and 1, 2 indicate more complex transitions.The bifurcation structure around 2 could not be fully detailed but findings are consistent with the picture provided in the inset.further, the shape and bubble pressure of the flat-tipped finger in the simulations approach those of pointed fingers, with pressure below the experimental data.The insets indicate that the overall finger shape upon which instabilities develop in the experiment is the pointed finger (blue) rather than the flat-tipped finger (red); see inset 4. As increases, it becomes progressively more difficult to distinguish these two modes; see inset 5.
In figure 10(b) we highlight the symmetry and stability of the steady solutions which were computed by imposing fixed * and to mimic the experiments.The different regimes are highlighted with shaded regions that are defined by the pushing/peeling transition, where the slope changes sign, and the viscous/elastic transition, estimated from 4(c).The solution structure is complex, consistent with that previously reported by [14], but does not always have a direct connection to the observed time-dependent behaviour.Consequently, we have not pursued a detailed study of this bifurcation diagram, other than to confirm that all bifurcations and changes of stability are consistent with standard theory.There are, however, connections that can be made between the bifurcation diagram and experimental observations.The round-tipped symmetric finger observed at low loses stability to a pair of emerging asymmetric fingers at a supercritical pitchfork bifurcation ( 1 ), which is consistent with the experimental transition to a flat-tipped asymmetric finger at = 0.1.The asymmetric fingers are superposed in this projection because they have the same bubble pressure.The unstable symmetric branch emanating from ( 1 ) is double-tipped, which provides a mechanism for intermittent tip-splitting near = 0.1 as seen in figure 8, panel 3.The finger may visit the vicinity of this unstable branch in the experiment intermittently due to background perturbations (figure 8(a), panel 3) or once following initial perturbation in the numerical system (figure 8(b), panel 3).
The pitchfork bifurcation 2 occurs at = 0.18 and the bifurcation scenario in its vicinity is consistent with the experimental transition from asymmetric flat-tipped fingers to pointed fingers, which accompanies the abrupt saturation of RW in §IV A. This signals elastic reopening, where * b is proportional to , and thus suggests that finger morphologies such as pointed and feathered modes are a feature of elastic vessel reopening [10,13,25].In the numerical model, the feathered modes emerge beyond the Hopf bifurcation 4 .Between 4 and a subsequent Hopf bifurcation, 5 , the pointed finger has one pair of unstable complex conjugate eigenvalues, while after 5 it has two pairs.However, in both regions, the frequencies of the eigenmodes do not match the typical oscillatory time scales of the feathered modes; see Figure 9.This suggests that these complex modes of propagation, discussed in §IV C 1, correspond to fully nonlinear dynamics.It is worth pointing that there are likely to be other solutions, in addition to those presented here.For instance, [14] have found multi-tipped Romero-Vanden-Broek [15,19,41,45] solutions, calculated originally for the rigid channel.These solutions were not stabilized by the presence of the elastic sheet.

V. CONCLUSION
We have explored the transition to complex pattern formation associated with two-phase flows in a Hele-Shaw channel, where the upper rigid boundary is replaced by an elastic sheet.This elastic upper boundary enables the collapse of the channel so that it adopts a prescribed non-uniform cross-sectional depth distribution and its cross-sectional area is reduced relative to the undeformed rectangular channel.Injection of air at constant flow rate into this initially-collapsed, liquid-filled channel leads to the propagation of an air finger, which reopens the channel by redistributing resident fluid ahead of its tip.We find experimentally that these propagating fingers exhibit unsteady dynamics for sufficient large level of initial channel collapse (small ) or large finger capillary number ( ), and that the characteristic feathered modes first identified by [10] are supported across a wide range of and .To capture these complex finger patterns with numerical simulations, we extend the depth-averaged model developed by [14] to include an artificial disjoining pressure term.This term is added to circumvent self-intersection of the air finger interface which prematurely terminates the simulations, and occurs when the applied local interface perturbation grows into a cleft that is sufficiently deep and narrow.The fact that self-intersection of the air-finger interface is not observed in the experiments suggests that its presence in simulations is likely a result of the approximations made when simplifying the three-dimensional liquid-film dynamics.
Although the inclusion of fluid film effects in the depth-averaged model is fundamental for both qualitative and quantitative agreement with experiments [14], we find that the results from numerical simulations are unaffected by our precise choice of film model, i.e., whether the film thickness is set at the finger tip and thus has a uniform thickness [38], or varies with the local height of the finger [14].This is because the choice of fluid film model only affects the system via mass conservation, so our results are unchanged provided we present them as a function of the capillary number rather than the imposed flow rate.
We started by quantifying the relative importance of viscous and elastic effects on the flow in the reopening region ahead of the finger tip.For this we compared the elastic system with the propagation of a finger into a rigid wedge [32,37].A rigid-wedge capillary number ( RW ) was estimated by converting the pressure gradient ahead of the finger tip, using the depth-averaged theory, into a dimensionless finger speed set by viscous and capillary forces.We used this metric to establish the fidelity of the simulations across the range of and studied, by demonstrating notable quantitative agreement between numerical and experimental values of RW .Comparing RW to , which is measured based on the finger speed, we find an abrupt saturation of RW at a threshold value of for sufficiently large initial collapse.This saturation of viscous dissipation signals a transition from reopening dominated by viscous and capillary forces to dominant elastic and capillary effects.We refer to these distinct reopening regimes as 'viscous' and 'elastic', respectively.We find that the exotic finger morphologies observed experimentally, such as pointed fingers and the feathered mode, are associated with an elastic regime of reopening.
A combination of experiments and numerical simulations was next used to explore the modes of finger propagation associated with the viscous and elastic regimes.In the viscous regime which corresponds to modest channel collapse and/or moderate values of , we recover the non-monotonic variation of bubble pressure (and width) with , which is characteristic of benchtop models of pulmonary airway reopening [16,30], and fundamentally distinguishes elastic-channel finger propagation from two-phase displacement flows in rigid Hele-Shaw channels [42].
The introduction of a disjoining pressure into the numerical model enables the study of the elastic regime, which occurs for increased channel collapse and/or .This is because the disjoining pressure prevents the self-intersection of the interface when the imposed perturbation evolves into a cleft, e.g., for ≥ 0.24 when = 0.53.Remarkably, the extended numerical model captures the destabilization of the finger into feathered modes of propagation, and this process helps to shed light on the rich variety of pattern formation observed experimentally.The feathered modes of propagation emerge following long asymmetric ( = 0.31) or symmetric ( = 0.32) oscillatory transients, which are associated with patterns reminiscent of the oscillatory fingers observed in the experiment for = 0.24 and 0.28, respectively.Furthermore, both experiments and simulations indicate that the feathered modes themselves may be transient.Increasingly fine indentation of the finger tip as it propagates indicates that the time-scale for tip-splitting at the interface is shorter than the time-scale for advection of the interface perturbation, and thus leads to ever refining fingering.These findings suggest that the unsteady patterns observed experimentally, characteristic of the elastic regime, may be continually evolving, with no evidence that these disordered modes converge to steady or periodic states.
Finally, we find excellent quantitative agreement between experiments and stable steady numerical solutions.Steady solutions also match the bounding envelopes of unsteady fingers obtained in the experiments and time simulations.This means that the steady bifurcation structure calculated numerically is sufficient to predict the bubble pressure.However, the steady bifurcation diagram does not inform unsteady finger propagation prevalent in the elastic regime, in that time-dependent modes which bifurcate from the steady solutions have different characteristic time scales from the feathered modes observed in the experiments and time simulations.Hence, in spite of the apparent simplicity of this system, the dynamics observed in the experiments and time simulations reflect complex nonlinear dynamics that require finite amplitude perturbations of the steady bifurcation structure to be attained.the closest agreement is obtained with a uniform film thickness distribution.Figure 13 shows that the steady finger solutions obtained using the two variations of the film model superimpose for the same capillary number.Steady results from both models, such as finger shape, bubble pressure and transition between modes of propagation, are in remarkable agreement when we use the capillary number as the control parameter.Hence, throughout this manuscript, we present all experimental and numerical results in terms of capillary number.However, since the capillary number is time-dependent, in the in the unsteady propagation scenario, the non-uniform film thickness model is the natural choice for time-dependent simulations.Therefore, all numerical results presented in this manuscript, apart from those in this appendix, are obtained with a choice of non-uniform film thickness.

FIG. 1 .
FIG. 1.(a) Schematic diagram of the experimental set-up.Adapted from [10].(b) Initial elastic-sheet profiles across the width * of the channel.Membrane height * , measured by the laser sheet depicted in (a), as a function of the scaled lateral coordinate * 2 / * is normalised by the channel depth * 0 .Labels refer to , a measure of initial collapse.(c) Top: image of air finger within the experimental region of interest (ROI).This image was captured by the top-view camera during the reopening of a channel with = 0.53 at flow rate * = 90 ml/min.Bottom: composite image produced by overlaying successive images from the experiment, taken at time intervals of 0.17 s.Time increases from left to right.(d) Instantaneous elastic-sheet height * mid-way between the channel walls ( * 2 = 0) as a function of axial coordinate * 1 during reopening [taken from the same experiment as (c)].The dashed line indicates the position of the air-oil interface at the tip of the air finger.

FIG. 2 .
FIG. 2. (a) Numerical domain in a frame of reference moving at the velocity of the finger tip.The rigid side walls are located at 2 = 0.5 and 2 = −0.5.The domain is truncated at 1 = − up upstream and at 1 = down downstream.The 1 coordinate of the finger tip is fixed at 0 and the finger width is .(b) Sketch of the transverse view of the channel along the line 2 = 2,tip that crosses the finger tip.The thickness of the fluid film is 1 ( ) and the thickness of the air finger is − 1 ( ) .The height of the elastic sheet at the finger tip is denoted by tip = ( 1,tip , 2,tip ).

FIG. 3 .
FIG. 3. (a) Example of a numerical time-dependent solution where the air-finger interface develops an indentation.One in every 25 nodes is shown with circles on the interface along with their normal vector.The magnitude of the vector ì corresponds to the distance between the points and .At the point , the interface-interface distance is equal to | ì |.(b) Example of an unsteady simulation where the finger develops an indentation in the absence of disjoining pressure.(c) The same unsteady simulation as in (b), and at the same time, but with the addition of disjoining pressure.

FIG. 4 .
FIG. 4. (a) Sketch of the fluid wedge in front of the finger.(b) Rigid-wedge capillary number (based on the pressure gradient in the fluid wedge ahead of the finger) as a function of the mean capillary number of the propagating finger, time-averaged during the propagation over the ROI.Circles indicate experimental data with error bars denoting standard deviations within the ROI and the black dashed line corresponds to the limit RW = .(c) Comparison between experimental data from (b) and numerical simulations shown with lines.The panels from left to right are for = 0.53, 0.60, 0.82 and 0.95, respectively.

FIG. 5 .
FIG. 5. Steadily propagating fingers for = 0.95.In experiments (1-5), the mean value of the capillary number over the ROI is = 0.01, 0.17, 0.45, 0.78, 1.23 respectively, while time intervals between the interfaces are 3.0, 0.3, 0.1, 0.05 and 0.05 s, respectively.The interfaces in red are the steady numerical solutions for the same capillary number as in the experiments.

FIG. 6 .
FIG. 6. Slightly collapsed channel, = 0.95.Plots of (a) the finger width and (b) the bubble pressure * b as functions of the mean capillary number .The solid black lines depict the numerical steady solutions, while the red symbols indicate mean experimental values with error bars corresponding to standard deviations within the ROI.In both figures, we identify a pitchfork bifurcation where the symmetric round-tipped branch bifurcates into an asymmetric flat-tipped branch.

FIG. 7 .
FIG. 7. Experimental images of fingers propagating at a fixed mean capillary number = 0.36 over the ROI.From top to bottom the initial level of collapsed increases.In experiments 1 to 5: = 0.78, 0.68, 0.56, 0.53, 0.43.The time interval between interfaces is the same for all experiments, 0.14 s.

FIG. 8 .
FIG. 8. (a) Time evolution of experiments performed at a fixed = 0.53 and constant flow rate that increases from 5 ml/min in experiment 1 up to 330 ml/min in experiment 10.The time interval between the interfaces are, from 1 to 10: 1.17, 0.67, 0.40, 0.17, 0.17, 0.17, 0.17, 0.10, 0.17, 0.13 and 0.10 s, respectively.The third to last interface in experiment 6 is missing due to a camera fault.See also Supplementary Movies.(b) Time evolution of our time-dependent simulations at a fixed = 0.53 and constant flow rate that increases from 3 ml/min in simulation 1 up to 160.3 ml/min in simulation 10.The time interval between the interfaces is 0.08 in the non-dimensional scale, which results in 1.486, 0.205, 0.134, 0.078, 0.062, 0.052, 0.043, 0.038, 0.037, 0.028 s from 1 to 10, respectively.The initial steady solution in blue and the perturbed profile at in red are the first two finger profiles in each time-sequence.See also Supplementary Movies.(c-d) Complete domain of the time-simulations from panel 8 and 9 are depicted in (c) and (d), respectively.
FIG. 8. (a) Time evolution of experiments performed at a fixed = 0.53 and constant flow rate that increases from 5 ml/min in experiment 1 up to 330 ml/min in experiment 10.The time interval between the interfaces are, from 1 to 10: 1.17, 0.67, 0.40, 0.17, 0.17, 0.17, 0.17, 0.10, 0.17, 0.13 and 0.10 s, respectively.The third to last interface in experiment 6 is missing due to a camera fault.See also Supplementary Movies.(b) Time evolution of our time-dependent simulations at a fixed = 0.53 and constant flow rate that increases from 3 ml/min in simulation 1 up to 160.3 ml/min in simulation 10.The time interval between the interfaces is 0.08 in the non-dimensional scale, which results in 1.486, 0.205, 0.134, 0.078, 0.062, 0.052, 0.043, 0.038, 0.037, 0.028 s from 1 to 10, respectively.The initial steady solution in blue and the perturbed profile at in red are the first two finger profiles in each time-sequence.See also Supplementary Movies.(c-d) Complete domain of the time-simulations from panel 8 and 9 are depicted in (c) and (d), respectively.

FIG. 9 .
FIG. 9. Time evolution of the 2 coordinate of a point on the interface at a fixed axial distance = 0.5 behind the finger tip.For (a) = 0.24, (b) = 0.28, (c) = 0.31 and (d) = 0.32.These values are extracted from the simulations depicted in figure 8(b).
Figure 10(a) compares bubble pressure * b in the experiments and steady numerical simulations as a function of the mean capillary number for = 0.53.As for weakly collapsed channels in §IV B, * b increases with indicating peeling fingers.

FIG. 10 .
FIG. 10.Plot of the bubble pressure * b as a function of the mean capillary number at the fixed level of collapse = 0.53.(a) Experimental fingers are plotted as symbols.Green circles represent round-tipped, red squares represent flat-tipped (symmetric or asymmetric), blue triangles represent pointed-tipped and purple diamonds represent feathered fingers.Filled symbols represent steadily propagating experimental fingers, while empty symbols represent unsteady ones.Error bars indicate standard deviations of the capillary number within the ROI.Steady numerical solutions are presented as solid lines following the same color code from the experiments.Insets are experimental images with steady numerical interfaces at the same superposed.They correspond to the numbered points in the graph.(b) Bifurcation diagram of the steady numerical solutions.Solid (dashed) lines represent stable (unstable) solutions.Black (blue) lines represent symmetric (asymmetric) solutions.The number of positive real eigenvaluesand complex eigenvalues with positive real parts , corresponding to instabilities, are indicated in brackets ( , ) for each solution branch.The relevant bifurcations are marked with red dots. 1 and 2 are pitchfork bifurcations, 1 -5 are Hopf bifurcations, 1 is a limit point and 1, 2 indicate more complex transitions.The bifurcation structure around 2 could not be fully detailed but findings are consistent with the picture provided in the inset.

FIG. 11 .
FIG. 11.Relationship between the transmural pressure and the level of collapse of the channel, .The red circles indicate the experimental data and the solid black line corresponds to the numerical simulation for (0) * 22 = 32 kPa.

Figures 14 and 15
Figures 14 and 15  show the results from experiments and steady simulations conducted at fixed levels of collapse = 0.82 and = 0.60, respectively.In each figure, part (a) shows the time evolution of the fingers for increasing values of and part (b) the bubble pressure as a function of .Finger propagation at = 0.82 is broadly similar to observations at = 0.95 described in §IV B, but here the fingers transition from round-tipped symmetric ( = 0.51) to flat-tipped asymmetric ( = 0.90).

FIG. 13 .
FIG. 13.Interfaces of steady solutions computed using the uniform (non-uniform) approach of the fluid film model plotted as the blue (black) solid line.The solutions are computed by fixing the value of capillary number: (a) = 0.06,(b) = 0.15,(c) = 0.23 and (d) = 40.
) is the local height of the channel and * = * up + * down = 25 * is the channel length.Throughout this manuscript, dimensional quantities are always starred, while dimensionless ones appear unstarred. *