Two-interface and thin filament approximation in Hele–Shaw channel flow

When a viscous fluid partially fills a Hele–Shaw channel, and is pushed by a pressure difference, the fluid interface is unstable due to the Saffman–Taylor instability. We consider the evolution of a fluid region of finite extent, bounded between two interfaces, in the limit the interfaces are close, that is, when the fluid region is a thin liquid filament separating two gases of different pressure. In this limit, we derive a second-order ‘thin filament’ model that describes the normal velocity of the filament centreline, and evolution of the filament thickness, as functions of the thickness, centreline curvature and their derivatives. We show that the second-order terms in this model, that include the effect of transverse flow along the filament, are necessary to regularise the instability. Numerical simulation of the thin filament model is shown to be in accordance with level-set computations of the complete two-interface model. Solutions ultimately evolve to form a bubble of rapidly increasing radius and decreasing thickness.


Introduction
For a standard model of Hele-Shaw flow in a rectilinear channel, consisting of a semi-infinite inviscid fluid region and a semi-infinite viscous fluid region separated by a single interface, the interface exhibits the Saffman-Taylor instability when the inviscid fluid displaces the viscous, and is stable when the viscous fluid displaces the inviscid (Saffman & Taylor 1958).In the absence of surface tension, exact solution methods exist that exhibit either finite time singularity or long-time finger formation (Howison 1986a,b); however, the problem is ill-posed.A regularisation such as surface tension is needed to stabilise sufficiently large wavenumber perturbations.In the presence of surface tension, solutions generically tend to a linearly stable travelling wave solution known as the Saffman-Taylor finger, although for very small surface tension, finite-amplitude noise in experiments and numerical simulations can lead to tip-splitting (Casademunt 2004).
The traditional Saffman-Taylor instability is normally studied on the assumption that the viscous fluid region extends infinitely far along the channel, so that there is only a single interface to consider.However, in reality the fluid region will only have finite extent, so that there are in fact two interfaces: one in which the viscous fluid is being displaced by the driving inviscid fluid, and one in which the viscous fluid is the one advancing (see Fig. 1).In this case, the force driving the fluid region is the pressure difference between the two inviscid fluids on either end of the fluid region.In the absence of surface tension, some classes of exact solutions to the two-interface problem have been found through use of special functions, in both channel geometries (Crowdy & Tanveer 2004;Feigenbaum et al. 2001) as well as annular geometries (Crowdy 2002;Dallaston & McCue 2012).The exact solutions in these studies, however, do not exhibit the interfaces becoming closely separated.Richardson (1982Richardson ( , 1996) ) similarly finds exact solutions for two-interface channel flow, some of which exhibit topological changes when two cusps form on different parts of the interface at the same point in space (none of these solutions correspond to the two semi-infinite regions of different pressure meeting, however).Farmer & Howison (2006) consider an approximate model where the two interfaces in an unbounded Hele-Shaw cell are very close, resulting in a thin filament of viscous liquid, and construct exact solutions to this approximate model (see Section 4.1).All of these zero-surface-tension models are ill-posed.
The two-interface Hele-Shaw model in the limit that the interfaces are close is of great importance as, even if the fluid region is not initially thin, the effect of the Saffman-Taylor instability will result in a thin fluid region or filament developing (see for example the experimental results in Ward & White (2011); Morrow et al. (2023), and Cardoso & Woods (1995), although the latter study focuses on a much less viscous thin annular region in a more viscous ambient fluid).This formation of a thin filament precedes the fluid 'bursting', at which point the two inviscid regions meet and the pressure rapidly equalises.The breakup of a thin viscous filament in a Hele-Shaw cell (with surface tension but in the absence of a driving pressure difference) has also been examined by the use of the lubrication approximation (Almgren 1996;Almgren et al. 1996;Constantin et al. 1993;Dupont et al. 1993;Goldstein et al. 1993Goldstein et al. , 1995Goldstein et al. , 1998)).The complicated, but self-similar break-up behaviour of the filament in particular (where the filament thickness goes to zero at a finite time and point in space) is detailed in Almgren et al. (1996).
In this article we consider two-interface Hele-Shaw flow in a channel, including the effects of both driving pressure difference and surface tension, with particular focus on the case in which the fluid interfaces become close together.In Section 2 we describe the full twointerface model and its stability.In Section 3 we derive a simplified model (the thin filament approximation) that applies when the two interfaces are close together, by applying the lubrication approximation (up to second order) in a coordinate system following the filament centreline.The application of lubrication theory in an evolving, curvlinear coordinate system is similar to that applied in Van De Fliert et al. (1995) and Howell (2003), although taking such an approximation to higher order is not standard.Our approximation also represents a generalisation of the model by Farmer & Howison (2006) by including both the effects of surface tension, as well as higher order terms, which include the effect of flow along the filament.In section 4 we demonstrate the importance of including these higher order terms by showing (both numerically and by constructing similarity solutions) that the leading order problem generically blows up in finite time, even in the presence of surface tension.We also find quasi-travelling wave solutions, which may be thought of as the analogue of Saffman-Taylor fingers, which generalise the 'grim reaper' exact solution found by Farmer & Howison (2006).In Section 5 we compute numerical solutions of the thin filament model including the higher order regularising terms.Solutions of the original two-interface model, found using a level set method, are seen to approach the thin filament model results as the resolution of the level set method is increased.We show the general behaviour of our thin filament model is not to tend toward a quasi-travelling wave solution, but instead develop a rapidly expanding 'bubble' of circular shape and decreasing thickness.

Formulation of the two-interface Hele-Shaw flow problem
2.1.Hele-Shaw flow equations We consider flow in a Hele-Shaw channel of nondimensional width 2 in the -direction.The fluid region Ω is bounded above and below by interfaces Ω U and Ω L , respectively.A nondimensional pressure difference  acts to push the fluid region in the positive -direction, while surface tension acts on both interfaces (see Fig. 1).The standard governing equations in nondimensional form are (Morrow et al. 2021, e.g.) with  the velocity potential,  n the normal velocity of each interface,  U and  L the curvatures of each interface,  is the nondimensional pressure difference, and  the nondimensional surface tension.In (2.1) the normal vector n of each interface is defined so that a flat interface has normal in the positive -direction, and the signs of curvatures  are chosen such that positive  implies a concave interface in the positive -direction (for example, both interfaces in Fig. 1 have negative curvature at  = 0 and positive curvature at  = ±).This choice of sign convention for both interfaces will simplify the thin filament derivation in the next section.
For definiteness, we consider no-flux boundary conditions on the channel walls: although many of our results, in particular the thin filament model derived in Section 3, do not rely on these conditions, and could equally apply to an infinitely long fluid region in an unbounded cell, or an annular fluid region represented by a closed curve (in that case, the constant pressure difference  would physically require injection or extraction of air into the interior of the cell).
The system (2.1) is derived from the dimensional system by introducing the length scale [], the pressure scale [ ], and defining the dimensionless velocity potential  in terms of the pressure : so that where   and   are the upper and lower inviscid fluid pressures,  is the plate separation, ℓ is the channel width,  the viscosity, γ is the surface tension, and [] is an arbitrary time scale.While one of the parameters  or  (if nonzero) could be scaled out by choosing an appropriate time scale, retaining both parameters aids in understanding the effects of the terms that arise in our later analysis.

Linear stability
The two-interface Hele-Shaw configuration has an exact base state where both interfaces are horizontal, and the fluid region moves upward at constant velocity  n = /ℎ 0 , where ℎ 0 is the distance between the two interfaces.To examine the stability of this configuration, write the system (2.1) in Cartesian (, ) coordinates, and define the upper and lower interfaces as  =  U (, ) and  =  L (, ), respectively.In addition to Laplace's equation for , the boundary conditions (2.1b)-(2.1d)are (2.3a) (for notational compactness we will use variable subscripts to refer to partial derivatives throughout this article; text or numerical subscripts, however, do not refer to partial derivatives).The base state is (up to arbitrary translation in ) represented by  L = (/ℎ 0 ) − ℎ 0 ,  U = (/ℎ 0 ), and  = (/ℎ 0 ) ŷ, where ŷ = ( − /ℎ 0 ) is the coordinate in the travelling frame.
To examine linear stability then, we impose perturbations on  L ,  U : On substitution into the boundary conditions, A perturbation with wavenumber  in the -direction takes the form Figure 2: Growth rate  of perturbations of given wavenumber  for the two interface Hele-Shaw model (2.5), the filament model (3.9) (second-order in the lubrication parameter), and the leading-order version of this model, which lacks terms corresponding to transverse flow (4.2).Each model has the parameter values  = 1, ℎ 0 = 0.2, and  = 0.1.All systems have two eigenvalues for each wavenumber.For the Hele-Shaw model and the filament model, the most unstable eigenvalue is stabilised at a cut-off wavenumber , and agree closely.The leading order filament model (4.2) is unstable for all wavenumbers, with the positive eigenvalue tending to a constant as  → ∞.
and , and the kinematic conditions result in an eigenvalue problem for : The eigenvalues of this system are thus The growth rate curves ( = ()) are shown in Fig. 2 (solid line) for parameter values  = 0.1, ℎ 0 = 0.2.One eigenvalue is negative for all wavenumbers , while the other is positive for a finite band of wavenumbers that ranges from zero up to a critical wavenumber.The presence of surface tension regularises the system by stabilizing the wavenumbers for large .We emphasise that the system is unstable no matter the direction of the pressure gradient (regardless of the sign of ), as each direction involves one of the interfaces moving in the unstable direction according to the Saffman-Taylor instability.
Our linear stability analysis here for a finite viscous fluid evolving in a Hele-Shaw channel is analogous to that presented for a radial geometry with two interfaces (Morrow et al. 2023).Generalisations and alterations to include different viscous fluids on either side of the interfaces have been conducted for both linear and weakly nonlinear frameworks (Gin & Daripa 2015;Anjos & Li 2020).
We close this section by noting some limiting behaviour of (2.5).For large ℎ 0 , in order to keep the speed of the base state  (1), we need to keep the driving pressure difference  =  (ℎ 0 ).In that case,  ∼ − 3 ±  /ℎ 0 as ℎ 0 → ∞.This limit agrees with the wellstudied single interface problem (an infinite body of viscous fluid), with the plus (minus) sign associated with the unstable (stable) direction of flow.Of a more particular interest here is the other limit of ℎ 0 ≪ 1.Again, supposing that  =  (ℎ 0 ) in order to keep the interface speed  (1), we have where for (2.7)  =  (1) means strictly of order one, while for (2.8)  =  (1/ℎ 0 ) means strictly of order 1/ℎ 0 .

Model in intrinsic coordinate system
In this section we derive an approximation of the two-interface flow by considering the thickness of the fluid region (that is, the distance between the two interfaces) to be small.This approximation will be valid when the separation of the two interfaces is smaller than the radius of curvature of either interface, but still larger than the plate separation (if the interface separation is the same order as the plate separation, the depth averaging that leads to two-dimensional Hele-Shaw flow (2.1) is no longer valid).We will start by converting the governing equations and boundary conditions into an intrinsic coordinate system (, ) ↦ → (, ) in which the filament is represented by a centreline curve x(, ) = ( x, ȳ), with  the arclength coordinate, and  the coordinate normal to this centreline; thus x = x + n (3.1) (e.g., Van De Fliert et al. 1995;Howell 2003).The coordinate system is represented schematically in Fig. 3.Here n = (− ȳ , x ) is the centreline normal (again we use the variable subscript  to represent the  partial derivative, for brevity).In this coordinate system we specify the upper and lower interfaces to occur at  = ±ℎ(, )/2, respectively, so that the normal thickness is ℎ.The upper interface is thus given by the curve and a (non-unit) normal to this interface (as opposed to the centreline) is then where s = ( x , ȳ ) is the centreline unit tangent vector, and  = x ȳ − ȳ x is the centreline curvature.
The velocity normal to the interface (not the centreline )  U is equal to the normal derivative Figure 3: A schematic of the coordinate system (3.1)used to derive the second-order lubrication model (3.8), which describe the velocity  n of the centreline and the evolution of the film thickness ℎ, in a direction normal to the centreline.In the derivation of the model, it is important to distinguish between the centreline normal n and normals to the interface (e.g., n U on the upper interface), and velocity components in these respective directions.
of the potential via the kinematic condition (2.1b): The component of the upper interface velocity in the centreline-normal direction,  nU , is then Similarly, on the lower interface, the velocity in the centreline-normal direction is From the interface velocities (3.2) we compute the centreline velocity and thickness evolution.Let  be a Lagrangian coordinate (such that  is constant at a point on the curve x moving in the direction normal to the centreline); then Here the notation / represents the Lagrangian derivative (holding  constant).Taking the centreline-normal component of each term in this equation, and noting that n/ is orthogonal to n, we find: with the result for the lower interface derived in very similar fashion as that for the upper.
Arranging for  n and ℎ/, and substituting (3.2), we then find expressions for the normal velocity and thickness evolution in the centreline-normal direction: where the notation )/2 represents a value averaged between the two interfaces.
As is standard in free-boundary lubrication problems, it is useful to be able to express the thickness evolution in a conservative form.To do so, we first note that Laplace's equation (2.1a), under the coordinate transformation (3.1), becomes 1 Integrating over  gives The first term on the right hand side of (3.5) represents the effect of dilation or compression of the film as its length changes, while the second term represents the contribution of flux along the filament.
The expressions (3.3) and (3.5) are exact, but we have not yet determined the potential .We now formally take the lubrication approximation by substituting  =  , ℎ =  , and  = , taking  ≪ 1.In this expansion we assume the curvature  =  (1) (or smaller), and the surface tension  =  (1).Unlike most standard lubrication models where the leading order behaviour in  is all that is required, in our case we will see that terms at order  2 are necessary to fully regularise the instability, and so we will expand to this order.
Writing  =  0 +   1 +  2  2 + . . .and substituting into Laplace's equation (3.4), we have Here  is again the centreline curvature, and  1 and  2 are order  and  2 corrections to the curvature on the upper interface, respectively: (the corrections on the lower interface are the same, with appropriate changes in sign).
Solving for each component of the potential  term by term gives where is the leading-order centreline velocity.Substitution of these expressions into the expansions of (3.3) and (3.5) results in expressions for the rescaled normal velocity  n =  −1  n and / = ℎ/: In (3.7b), (3.7c), the order of terms in the lubrication parameter  is explicit.We can of course formally remove the scaling by returning to ℎ,  variables, resulting in with  0 = ( + 2)/ℎ the leading-order velocity, and  1 =  −1  1 ,  2 =  −2  2 the rescaled versions of (3.6): (3.8c) The higher order terms have been included in the above derivation as they are needed to fully regularise the instability.This will become apparent when we examine the behaviour of the leading order system in Section 4. In particular, the higher spatial derivative of thickness ℎ resulting from the interface curvature correction term  1 plays a crucial role.In the case where driving pressure  and centreline curvature  both vanish, (3.8) reduces to the well-known thin film equation , where ℎ = ℎ(, ), and arclength  becomes time-independent.This equation has been studied extensively in the context of droplet breakup in Hele-Shaw flow (Almgren 1996;Almgren et al. 1996;Constantin et al. 1993;Dupont et al. 1993;Goldstein et al. 1993Goldstein et al. , 1995Goldstein et al. , 1998)).In our case, this fourth order spatial derivative term stabilises high wavenumber perturbations, as is shown in the following stability analysis of the filament model.

Stability
As with the full two-interface model (2.1), the thin filament approximation (3.8) has an exact solution comprising a straight filament of uniform thickness ℎ 0 moving upward with speed /ℎ 0 .To test the stability of this straight filament, we write the centreline x(, ) in Cartesian coordinates as  =  (, ), write thickness ℎ as a function of  and , and perturb the straight filament: In the linear approximation, / = / +  (  2  ),  =  +  (  2  ) and  =    +  (  2  ).On substituting these expansions into (3.8)we obtain the eigenvalue problem The eigenvalues are thus (3.9) We note that the linear stability of the filament model (3.8) does not agree exactly with the full model, even for zero surface tension.When  = 0, the error in the filament model is due to the final term  2 ℎ 2 0  6 /144, which arises from the product of the two higher order correction terms in the off-diagonal terms in the determinant.This term is (implicitly) of order  4 , which is of higher order than the model is accurate to.If further terms were taken in the lubrication expansion, they would contribute further order  4 terms which would presumably cancel with the term present here.
For nonzero surface tension , (3.9) has the same limiting behaviours (2.6) and (2.7) as the full model, but differs slightly with (2.8).In the latter case, both (2.5) and (3.9) give  =  (1/ℎ 3 0 ) as ℎ 0 → 0 for  =  (1/ℎ 0 ), but with a different prefactor.Thus, our thin filament approximation recovers the same leading-order linear stability behaviour as the full model for small and moderate wave numbers, which is all we can expect from such a lubrication model.For very large wave numbers (i.e., very small wavelengths of perturbation), with  ≫ 1/ℎ 0 , while the scalings between the eigenvalues of the full and filament model may differ, these modes of perturbation decay very quickly and therefore any differences in the models are of no practical consequence.
In Fig. 2 we compare the eigenvalues of the thin filament approximation (3.9) against those of the full problem (2.5).For even moderately small thickness (ℎ 0 = 0.2) and small surface tension ( = 0.1), the agreement is excellent in the region in which eigenvalues are positive; the difference for negative values occurs for much larger  than depicted.These linear stability results do indicate, however, that the approximation will not be valid if surface tension is zero (or much smaller than the filament thickness), which is unsurprising given the implicit assumption that  =  (1) in the lubrication parameter .

Properties of the leading-order model
In this section we consider the properties of the leading-order version of (3.8), that is, neglecting terms of order  2 : where  is the centreline curvature.In the leading-order model (4.1), there is no fluid flow tangent to the centreline; the change in filament thickness is due purely to the stretching or compression of the filament.
Potential issues with the leading-order model can first be observed in the linear stability analysis of a flat filament.The stability analysis of the leading-order model (4.1) is the same as for the higher order model with the higher-order terms absent; the eigenvalues of this leading-order system are thus For small , (4.2) coincides with the leading-order behaviour (2.6), while for  =  (1), (4.2) is no longer a reasonable approximation for the full model.Indeed, unlike in (3.9),where both eigenvalues becomes negative for sufficiently large wavenumber , in (4.2), one eigenvalue tends to a positive, constant value as  → ∞: (see Fig. 2).Thus the system is (significantly) unstable to perturbations at arbitrarily small spatial scales, even in the presence of surface tension.While this analysis does not suggest the leading-order problem is technically ill-posed if  > 0 (the eigenvalues do not become arbitrarily large), this large- behaviour strongly suggests that singularities in curvature will generically form (see for example Dallaston & McCue 2014).We will confirm the self-similar formation of curvature singularities in Section 4.2.

The unregularised (zero surface tension) leading-order model
In the case  = 0, the leading order model (4.1) reduces to that considered by Farmer & Howison (2006).In this case the problem is indeed ill-posed; the eigenvalues (4.2) are  = ±(/ℎ 0 ), as are the eigenvalues of the full two-interface problem (2.5), so one eigenvalue is arbitrarily large as  → ∞.Solutions that exhibit this ill-posedness may be constructed by showing that the centreline is given by level curves of a harmonic function.We summarise the approach of Farmer & Howison (2006) and further examine this approach here.
Assuming the filament centreline ((, ), (, )) and the thickness ℎ(, ) are parametrised by a Lagrangian coordinate , then For  = 0, the evolution of ℎ (4.1) is equivalent to conservation of mass between a point  and reference point  0 .That is, we may define an area function (), such that is constant in time on a point on the centreline moving with its normal velocity.Choosing to scale time such that  = 1, we arrive at the following: These are Cauchy-Riemann equations relating (, ) to ( , ).This line of argument was used by Farmer & Howison (2006) to demonstrate that  = +i must be an analytic function of the complex spatial variable  =  + i; thus, for a given time , the centreline is the level curve of the harmonic function  (, ).
Given the definition of  (4.3), the thickness ℎ may be calculated by: This thickness will go to zero at a critical point   where  ′ (  ) = 0.As this is also a point where the conformal map between  and  ceases to be smooth, we expect to see a singularity in the curvature in the centreline there.The preimage of the straight line  =  + i  that passes through the critical point is the centreline in the -plane; assuming  ′′ (  ) ≠ 0 the centreline must therefore have a corner with an angle of /2 at   (and not a cusp, as suggested by Farmer & Howison (2006)).If the initial condition is such that  ′′ (  ) also vanishes but the third derivative is nonzero, the corner angle is /3, and so on.
As an example, consider an initial condition with the centreline on  = 0 and initial thickness given by ℎ(, 0) = [1 −  cos ], with  > 0 and 0 <  < 1.This initial condition corresponds to an initially horizontal filament that is thinner near  = 0. We thus have  = [ −  sin ] at  = 0 (determined by choosing our reference point  0 to lie on  = 0), and, analytically continuing into the complex plane: Taking the imaginary part, we find that the centreline location is given implicitly by The critical point occurs for   = i cosh −1 (1/) and time The centreline profiles of this solution, along with the upper and lower interfaces (found by adding and subtracting half the thickness ℎ (4.5), respectively, in the normal direction) are plotted in Fig. 4a, showing the formation of the /2 angle as  →  −  .As an example of an initial condition that leads to a non-/2 angle, consider Again, the centreline is determined by taking the contours of the imaginary part of this function.In this case, since  ′ (  ) =  ′′ (  ) = 0, the corner angle at the critical point where ℎ → 0 is /3.Profiles of this solution are shown in Fig. 4b. Figure 4: Exact solutions to the thin filament equation in the absence of surface tension.(a) the solution (4.7) that evolves from an initial condition that results in corner formation with the generic angle of /2, and (b) the solution (4.8) that evolves from an initial condition that results in corner formation with a non-generic angle of /3.The parameter values are  = 0.1,  = 0.2.In both of these examples, the centrelines are plotted in solid blue while the black dashed lines are the upper and lower interfaces (found by adding or subtracting half the thickness ℎ in the normal direction).(c), the profiles, and (d) the thickness, for the example (4.9) with  = 0.2 and  = 0.5, which forms a zero angle cusp with infinite thickness.
The previous two exact solutions show only a subset of the possible singular behaviours of the ill-posed zero-surface tension model; as centrelines are given by level curves of a harmonic function  = ℑ(()), with thickness given by (4.5), a variety of other solution behaviour is possible.One such possibility is a zero-angle cusp that results from a square root singularity, an example of which is the complex function This function has a square root singularity at  = 0, as well as the appropriate periodicity.The level curve that passes through this singularity occurs at  = 0, so as  → 0 − , a zero angle, backward-facing cusp forms on the interface.The thickness ℎ, which is given by | ′ ()|, becomes infinite at the cusp, so that the interface at the cusp becomes stationary.The profiles and thickness corresponding to (4.9) are shown in Fig. 4c-d, respectively, for the parameter values of  = 0.2 and  = 0.5.

Singularity formation in the leading-order model with surface tension
We now further consider the leading order-model (4.1) in the presence of surface tension.As eigenvalues (4.2) are positive for arbitrarily large wavenumber, we expect the possibility of singularities in curvature, and this does indeed generically occur, as we demonstrate through both the construction of similarity solutions, and corroborate with numerical simulations.
To establish the existence of curvature singularities we perform a self-similar analysis.This analysis is simplest to carry out in Cartesian coordinates.Let the thickness ℎ = ℎ(, ) be a function of  and , and define the centreline to be given by the function  =  (, ); we then have where / is the time derivative with  held constant.Using (4.1) then results in Further simplification is achieved by defining the filament thickness in the  direction h(, ) = ℎ √︁ 1 +  2  (as opposed to the thickness ℎ in the centreline-normal direction), which results in the system of equations . (4.10)By introducing h, (4.10) has been written in conservative form, which highlights the fact that mass is indeed conserved in this system.
Assume a singularity occurs at a time  =   at  =   , and let: where the similarity exponents , ,  are to be determined.Furthermore, we assume the profile and thickness are symmetric in .Assuming that  > 1, the dominant term in the velocity   is  0 =  0 (  ).Thus, on balancing terms we find  = −1 and  = 2 − 1, with  being undetermined (a second-kind self-similarity).Given  >  (which we will check for consistency after the fact), the dominant terms in (4.10) become These equations can be further scaled to remove  0 and .
Let  = F′ and eliminate F, so that  satisfies the second-order equation For symmetry we require  to be odd in ; thus, when  = 0,  = 0, which is therefore a singular point of the ODE.By expanding near this point we find that the similarity exponent  is uniquely specified by the first odd power in the expansion greater than unity: where  is an arbitrary scaling factor.Expanding near the singular point results in  = /( − 1) > 1, which is consistent with the assumptions made on the similarity exponents above.
Given  = /( − 1), the equation for  can be solved implicitly by letting  be the independent variable, which allows us to construct parametric solutions for F and Ĝ: While  can be any odd number ⩾ 3, dependent on the initial condition, the most generic case will be  = 3, in which case  = 3/2 and  = 2.The constant  is arbitrary, due to scale invariance in the equations (4.10) in the limit that curvature dominates over the driving pressure .In a given case, the scale , velocity  0 , and exponent  will all depend on the initial condition of the time-dependent problem.
We provide numerical evidence for the curvature singularity formation by numerically solving the system (4.10).This computation is performed in MATLAB using finite difference discretisation along with MATLAB's ode15s algorithm for time-stepping.Parameters  = 1,  = 0.5 and an initial condition of  (, 0) = − cos() and ℎ(, 0) = 1 is chosen in order to start with high curvature near  = 0, where the singularity will occur.
In Fig. 5 we plot the results of this numerical computation.The singularity time   is estimated by fitting a straight line through ℎ max () −1 = (max  ℎ(, )) −1 , which occurs at  = 0, and should (according to the analysis above) go to zero linearly.The centreline velocity  0 () at the maximum thickness is observed to tend to a nonzero constant, from which  0 is estimated.Scaling the profiles near the singularity time into similarity variables , F, Ĥ, we observe collapse.The exact similarity profiles (4.12), with a suitable fitted value of the arbitrary constant , match well with the numerical profiles.
The curvature singularity exhibited by this model is weaker than a corner singularity, in that as  > , the first derivative goes to zero in the neighbourhood of the singularity, even while the curvature becomes unbounded.It is also interesting to note that the singularity is not dependent on (and does not require) the driving pressure ; its cause is the presence of surface tension pulling regions of high positive curvature inward, concentrating the thickness at a single point.It is thus only present as the leading-order model (4.1) has a surface tensiondependent velocity, but no regularising term that penalise higher derivatives of thickness, as appears in the second-order model in (3.8b).

Quasi-travelling wave solutions
The system (4.10), while insufficiently regularised to simulate the full dynamics of a thin filament, does exhibit a quasi-travelling wave solution of the form (in Cartesian variables) where  0 is a finite time.A solution of this form is similar to, but not exactly, a travelling wave, as the centreline has a fixed shape but moves to infinity (with speed unbounded) as  →  0 , while the thickness linearly decreases to zero.The parameter  is the analogue of a travelling wave speed.Solutions to (4.10) of this form generalise the 'grim-reaper' solutions to the zero surface tension problem found in Farmer & Howison (2006), and are also candidates for asymptotically valid solutions to the higher-order system (3.8),since as the thickness goes to zero, we would expect the higher-order terms to vanish at faster rate than the leading-order terms.We will see in Section 5 that these quasi-travelling waves do not appear to be attractors, but will compute them here for completeness.On substitution of the ansatz (4.13) into the leading-order model in Cartesian form (4.10), and scaling the variables according to The singularity time (large dot) is found by fitting a straight line approximation to the reciprocal of the maximum thickness near the singularity time, while (f) the speed of the centreline  0 =   (0,   ) (large dot) at the singularity is similarly estimated.The five smaller points in (e,f) are the times at which the scaled profiles are plotted in (c,d).
we obtain This equation may be solved numerically directly, but to further simplify we cast it into the following curvature-angle formulation.Let  = − tan −1 F′ be the angle between the tangent to the centreline and the -axis (counting positive for negative  ′ ), and scaled curvature  =  −1  = F′′ /(1 + F′2 ) 3/2 (see Fig. 6a).We then obtain the first-order equation d d For a semi-infinite curve, the appropriate interval is −/2 <  < /2, where the nose at  = 0 is a singular point, at which we require  (0) = −1.In the limit of the tail of the travelling wave  → /2, we must have  ∼  − /2 → 0. We solve equation (4.15) numerically for different values of σ, and reconstruct the x, F coordinates, using The σ = 0 and σ → ∞ limits are amenable to exact solutions.If σ = 0 then  =  0 = − cos(), which is the 'grim reaper' solution found by Farmer & Howison (2006) (also relevant for the zero-surface-tension Saffman-Taylor finger (Saffman & Taylor 1958) and for curve shortening flow (Angenent 1991)), equivalent to a centreline profile F = ln(cos x).On the other hand, as σ → ∞ the equation for  becomes linear, and has exact solution Plots of these quasi-travelling wave solutions are shown in Fig. 6b.The effect of σ, and thus of surface tension, is weak, as all solutions are bounded between the σ = 0 and σ → ∞ limits.The scaled curvature at the nose is required to be −1 in all cases, and thus returning to the unscaled system, we expect the curvature at the nose to be −1/, that is, inversely proportional to the wave speed parameter.
The corresponding scaled thickness in the  direction is H = sec 2  (1 + σ).The timedependent thickness in the centreline-normal direction is thus At a fixed time  0 then, the thickness is unbounded in the tails as  → ±/2.If, however, we consider a fixed value of  (in the non-travelling coordinate system) as the singular time  0 is approached, then from (4.13) and (4.16), the corresponding value of  behaves as Thus ℎ =  (e −/ ) is bounded for a fixed value of  as  →  0 .The exponential decay in  of the thickness also means that a quasi-travelling wave solution does not require infinite mass to form, even as it length is unbounded in the -direction as  →  0 .
In Section 5 we will simulate the second order thin filament model, including with initial conditions close to a quasi-travelling wave solution.These simulations will demonstrate that these quasi-travelling wave solutions do not appear to be stable.

Numerical simulation and comparison
In this section we describe numerical simulations of the second-order thin filament model (3.8), and compare against simulations of the full two-interface problem (2.1).
Our solution method for the thin filament model (3.8) is a front tracking method, whereby the centreline is represented by  points x  = (  ,   ),  = 1, . . ., .The thickness also has a value ℎ  at each point.At a given time, the normal vector, curvature, and arclength derivatives of ℎ and related quantities in (3.8) are calculated using a central finite difference scheme.The points are then moved in time in the normal direction with velocity given by (3.8a) (consistent with (3.8b) being Lagrangian evolution equation for ℎ) using MATLAB's ode15s.Since moving points with normal velocity results in highly unevenly spaced points on the centreline, we remesh onto an evenly spaced grid when the ratio between minimum and maximum node spacing drops below a threshold value.We typically use 2 000-5 000 points (increasing this resolution does not have an appreciable effect in the simulations reported here) and remesh when the minimum-to-maximum node spacing is less than 0.8.

Validation against solution of the two-interface problem
To validate solutions of the thin filament model (3.8) it is valuable to compare it to solutions of the full two-interface system (2.1), which is a more challenging numerical problem.We solve (2.1) using the numerical framework proposed by Morrow et al. (2021Morrow et al. ( , 2023)), which we briefly summarise here.The framework is based on the level set method, in which we represent each interface,  L and  U , as the zero level set of the associated level set functions  L and  U .Each level set function is chosen such that the viscous fluid will occupy the region where both  L and  U > 0; otherwise the region is filled with inviscid fluid.Both level set functions are updated according to where and n L = ∇ L /|∇ L | and n U = ∇ U /|∇ U | are the unit (outward) normals.We discretise the spatial derivatives in (5.1) using a second-order essentially non-oscillatory scheme and integrate in time using second-order total-variation-diminishing Runge-Kutta with time step Δ = 10 −5 .We perform simulations on the computational domain − ⩽  ⩽  and −0.5 ⩽  ⩽ 4, which is discretised into equally spaced nodes with mesh size Δ = 2/750.Simulations are concluded when the minimum distance between the two interfaces is less than 3Δ.Further, we periodically perform reinitialisation to maintain both  L and  U as signed distance functions.The details of this reinitialisation procedure are described in Morrow et al. (2021).We solve (2.1a) for the velocity potential  via a finite difference stencil.Following Chen et al. (1997), we modify the stencil at nodes adjacent to either interface, corresponding to where  L or  U changes sign, by imposing a ghost node on the interface to incorporate the appropriate dynamic boundary conditions (2.1c) and (2.1d).Here,  L = ∇ • n L and  U = ∇ • n U .By solving for , we can compute  L and  U from (5.2).These choices of  L and  U satisfy the kinematic boundary conditions (2.1b) and gives a continuous expression for  L and  U in the region occupied by the viscous fluid x ∈ Ω.However, to solve (5.1), we require expressions for  L and  U over the entire computational domain.To extend our expressions for  L and  U into the gas regions, we follow Moroney et al. (2017) by solving the biharmonic equation (5.3) By doing so, we obtain smooth continuous normal velocities over the entire computational domain, allowing us to solve (5.1) for  L and  U .
In Fig. 7 we compare the results of the filament model and full two-interface model, with initial conditions:  U (, 0) = 1/24 − 0.00375 sech 2 ,  L (, 0) = −1/24 + 0.00375 sech 2 . (5.4) These initial conditions correspond to an initial centreline location of  = 0 and an initial thickness ℎ(, 0) = 1/12−0.0075sech 2 , which is an almost flat filament with a small thinner region near the centre of the channel,  = 0.As demonstrated in Fig. 7a, the agreement between the two methods is initially very good, and only starts to disagree quantitatively when the filament thins near the central region, leading to a large increase in velocity.This is to be expected, as the level set method will become inaccurate when the filament thickness becomes too thin to capture accurately using a level set function on discretised mesh.To demonstrate this point, in Fig. 7b we plot the minimum thickness ℎ min () against time, for level set computations of increasing resolution.We see that the level set result does approach that of the filament model as the resolution is increased (of course, we do not expect the results to be identical due to the approximation made in deriving the thin filament model (3.8), but this error is much smaller than the achievable discretisation error in the level set simulation).Ultimately in the level set simulations, the thickness saturates at a value dependent on the resolution.This limitation means the level set simulations ultimately underestimate the speed at which the filament advances when it becomes very thin, as can be observed in the profiles shown in Fig. 7a.

Numerical results for late times
Here we use the front-tracking scheme to solve the thin filament model for later times, in the regime where the thickness becomes too small for the level set method to accurately resolve.We use the initial condition which is the same as the initial condition of the exact solution (4.7) depicted in Fig. 4b.In order to observe the effect of different surface tensions , we run simulations for  = 0.1 (Fig. 8), as well as  = 0.05 and  = 0.2 (Fig. 9).7: (a) A comparison between the interfaces predicted from numerical solution of the filament model, and those from numerical solution of the full problem using the level set method described in Section 5.1.The initial condition is given by (5.4) with pressure and surface tension parameters  = 1 and  = 0.1, respectively.For clarity, only the interfaces of the filament model, found from the actual variables of centreline and thickness, are plotted.(b) A comparison of the minimum thickness ℎ min (which occurs at the nose  = 0 between the filament model and level set simulations, as the resolution of the level set simulation is increased.For smaller Δ, the level set method tends to the thin filament result, while the filament thickness in the level set method saturates at a value that depends on the numerical resolution.The points marked on the filament curve correspond to the times at which the profiles are plotted in (a).
In the  = 0.1 simulation (Fig. 8), the filament bulges outward in the centre near  = 0, where the thin filament is initially thinner (and so the filament moves faster).This bulge becomes a 'bubble' that rapidly expands in radius while the thickness rapidly decreases.The majority of the fluid is pushed into the outer regions of the channel.At a finite time, the bubble intersects the channel walls at  = ±.Unlike a fluid-filled Hele-Shaw cell, there is nothing to prevent the filament reaching the channel walls; this does not correspond to singular behaviour in the mathematical model, but physically represents an area of gas of lower pressure being trapped by the filament.To further understand this behaviour, we note that (3.8) has as a solution a perfectly circular bubble of radius (), that evolves according to (5.5) Here  is a constant that depends on the initial thickness.In such a solution, the radius (if initially greater than 2/) grows exponentially, but does not exhibit a finite time singularity.It may be that this is the behaviour that is governing late stages of evolution of the filament depicted in this section.In Fig. 8e we plot the curvature  nose at the nose, or front, of the bubble (where  = 0), while in Fig. 8f we plot / nose .The curvature initially grows in magnitude (in our convention the curvature is negative in the bubble) when the bulge initially grows, but then rapidly decreases in magnitude at the time when the bubble expands outwards.When this happens, the curvature tends to a uniform state in the bubble region (/ nose → 1), implying the convergence to a circular shape.While we have only plotted the interfaces up to the time at which it intersects the channel wall, mathematically the simulation continues for longer, and the bubble becomes more circular in shape.Ultimately, the thickness becomes so small and the velocity so large that the numerical method no longer converges.This behaviour is also seen for smaller and larger surface tension values  = 0.05 and 0.2, as depicted in Fig. 9.The notable effect of different surface tension is that the bulging region that forms the 'neck' of the bubble is smaller or larger in width, as the surface tension is smaller or larger.

Initial condition near a quasi-travelling wave
In addition to the previous near-flat initial conditions, we demonstrate the instability of the quasi-travelling wave solutions considered in Section 4.3 by choosing an initial condition close to one such solution.For simplicity we use the 'grim reaper' exact solution for zero surface tension  0 = − cos , and a scale factor  = 1, for which  = log(cos ), ℎ =  cos  (the exact solution is a good approximation for small , and although we do not show the results here, we have observed the same behaviour from numerically constructed initial conditions by solving (4.15) for σ > 0).In terms of the quasi-travelling wave solution, the factor  in the thickness ℎ is arbitrary as it corresponds to translation in time.As the travelling wave is semi-infinite, we must of course approximate it by a finite curve that respects the required boundary conditions at  = ±.We thus choose the closely related 'hairclip' curve, defined by (5.6) The larger the parameter  is made, the more elongated the initial finger becomes.We note that the thickness will become large on the sides of the initial finger, representing the blow-up in thickness of the quasi-travelling wave in the limit  → /2 discussed in Section 4.3.We depict the result of such an initial condition in Fig. 10, for the values  = 0.05,  = 20,  = 0.1.The regions with large thickness on the sides of the finger correspond to the exponential growth in thickness required for the travelling wave solution as  → −∞ (see Section 4.3); the filament in this region does not evolve to a great extent over the simulation.At the front, while the finger does initially propagate forwards, the nose bulges outward and  becomes more circular in shape, leading to the same late-time behaviour as seen for the more general initial conditions.Figure 10: The evolution of a filament that starts near the ( = 0 limit) of the quasi-travelling wave described in Section 4.3 (the precise initial condition is (5.6)).(a) the quasi-travelling wave appears unstable, with the front evolving toward the expanding circle seen for general initial conditions.(b) The curvature at the nose of the filament ( = 0) starts at −1 but decreases over time as the filament expands.The dots marked in (b) correspond to the profiles plotted in (a).

Discussion
In this paper we have developed a simplified but highly accurate second-order lubrication flow model (3.8) that describes two-interface Hele-Shaw flow very well in regions where the thickness of the fluid region becomes small.Due to the instability of one of the interfaces, this limit is one that is generally reached, even if initially the thickness is not very small.By examining the generic singular behaviour of the leading-order model, even in the presence of surface tension (Section 4), we have shown why a second order-model is necessary to represent the dynamics of the full problem.In particular, the fourth order spatial derivative term arising from the difference in curvature between the upper and lower interfaces is necessary to regularise the system.Although we have included all terms formally of order  2 in our model and simulations, we have also observed that the leading order model along with the addition of only this regularising term (proportional to [ℎℎ  ]  ) will behave in a qualitatively similar manner (with small quantitative differences).We have not shown these results here for the sake of brevity.
Here we note some clear differences between the instability of a thin filament, and the classical Saffman-Taylor instability in a semi-infinite fluid region that results in the Saffman-Taylor finger with (in the small surface tension limit) width half that of the channel.One difference is that the thin filament model does not feel the effects of the channel wall away from the thick neck regions, and thus the rapidly expanding bubble may intersect the channel walls at finite time with no breakdown in the mathematical model.Physically speaking, this phenomenon would correspond to trapping a part of the lower-pressure gas inside the fluid region.In addition, as the walls have no strong effect, the orientation of the filament motion to be mainly in the positive -direction is a somewhat artificial consequence of the initial condition.For this reason it may be more natural to consider the fluid in an unbounded Hele-Shaw cell, with the length scale set by the initial thickness.
The specifics of the late stages of evolution of the filament (depicted in Section 5), wherein the thickness becomes very small and the velocity correspondingly large, is not resolved.In order to understand the dynamics at late times of this system, an analysis of the stability of the quasi-travelling wave solutions, and the stability of the expanding circular bubble (5.5) to non-radially symmetric perturbations, would be very valuable.We observe that in our model, the filament does not appear to exhibit finite-time 'bursting' behaviour, that is, the thickness does not go to zero at an isolated finite point in space and time.In the case of self-similar breakup in the manner of the unforced lubrication equation described in Almgren et al. (1996), the thickness goes to zero at a point where the curvature becomes infinite, while in our solutions depicted in Section 5, the thickness becomes small in the expanding circular region where the curvature is also decreasing, while the curvature is largest in magnitude in the neck regions, where the thickness does not decrease.Bursting in physical systems is likely to require fully three dimensional effects to explain (that is, when the filament thickness becomes the same order as the separation between plates in the Hele-Shaw apparatus).Once these two length scales are comparable, the filament model, and indeed the two dimensional Hele-Shaw model, are no longer valid.The dynamics of the filament model studied in this work indicate that such a regime will generically be approached very rapidly in the form of the expanding bubble depicted in Section 5.

Figure 1 :
Figure 1: A finite fluid region in a Hele-Shaw channel.The fluid is driven in the positive  direction by a pressure difference.

Figure 5 :
Figure 5: Curvature singularity formation in the leading-order Hele-Shaw filament model (4.10) for initial condition  (, 0) = − cos(), and h(, 0) = 1.(a) centreline profiles  (, ) approaching a curvature singularity and (b) filament thickness becoming unbounded at singularity time of   ≈ 0.77.These profiles (solid lines) collapse onto similarity profiles (c) F () and (d) Ĥ (), which asymptotically match the exact similarity solution (4.12) (dotted lines) as  →   .The scaling factor  ≈ 0.3 is fit to the profiles.(e)The singularity time (large dot) is found by fitting a straight line approximation to the reciprocal of the maximum thickness near the singularity time, while (f) the speed of the centreline  0 =   (0,   ) (large dot) at the singularity is similarly estimated.The five smaller points in (e,f) are the times at which the scaled profiles are plotted in (c,d).

Figure
Figure7: (a) A comparison between the interfaces predicted from numerical solution of the filament model, and those from numerical solution of the full problem using the level set method described in Section 5.1.The initial condition is given by (5.4) with pressure and surface tension parameters  = 1 and  = 0.1, respectively.For clarity, only the interfaces of the filament model, found from the actual variables of centreline and thickness, are plotted.(b) A comparison of the minimum thickness ℎ min (which occurs at the nose  = 0 between the filament model and level set simulations, as the resolution of the level set simulation is increased.For smaller Δ, the level set method tends to the thin filament result, while the filament thickness in the level set method saturates at a value that depends on the numerical resolution.The points marked on the filament curve correspond to the times at which the profiles are plotted in (a).

Figure 8 :
Figure 8: Evolution of a filament with initial position  = 0 and thickness ℎ = 0.2(1 − 0.1 cos()), with  = 1 and surface tension  = 0.1.(a) the centreline profiles over time (solid) and the thickness (dashed lines), and (b) the thickness against , show the initially thinner part of the filament bulge outward into a 'bubble', while the bulk of the fluid is driven out to the edges of the filament.Ultimately the profile intersects with the channel boundary at  = ±.The (c) maximum normal velocity  n max () and (d) minimum thickness ℎ min () appear to become unbounded and go to zero, respectively.(e) The curvature at the nose ( = 0) over time, initially increases in magnitude, then rapidly heads toward zero as the bubble expands.(f) the curvature scaled by the curvature at the nose, showing the curvature tending to a constant over the bubble.The dots marked in (c, d, e) correspond to the times at which profiles are plotted in (a, b, f), and the arrows in (a, b, f) indicate the direction of increasing time.

Figure 9 :
Figure 9: Evolution of a filament with initial position  = 0 and thickness ℎ = 0.2(1 − 0.1 cos()), for (a,b)  = 0.05,and (c,d)  = 0.2.The formation of a circular bubble occurs in a similar fashion as happens for  = 0.1, with the major effect being the width of the neck region that precedes the circular bubble.