Surface-tension-driven evolution of a viscoplastic liquid coating the interior of a cylindrical tube

One mechanism for airway closure in the lung is the surface-tension-driven instability of the mucus layer which lines the airway wall. We study the instability of an axisymmetric layer of viscoplastic Bingham liquid coating the interior of a rigid tube, which is a simple model for an airway that takes into account the yield stress of mucus. An evolution equation for the thickness of the liquid layer is derived using long-wave theory, from which we also derive a simpler thin-film evolution equation. In the thin-film case, we show that two branches of marginally-yielded static solutions of the evolution equation can be used to both predict the size of initial perturbation required to trigger instability and quantify how increasing the capillary Bingham number (a parameter measuring yield stress relative to surface tension) reduces the final deformation of the layer. Using numerical solutions of the long-wave evolution equation, we quantify how the critical layer thickness required to form a liquid plug in the tube increases as the capillary Bingham number is increased. We discuss the significance of these findings for modelling airway closure in obstructive conditions such as cystic fibrosis, where the mucus layer is often thicker and has a higher yield stress.

One mechanism for airway closure in the lung is the surface-tension-driven instability of the mucus layer which lines the airway wall. We study the instability of an axisymmetric layer of viscoplastic Bingham liquid coating the interior of a rigid tube, which is a simple model for an airway that takes into account the yield stress of mucus. An evolution equation for the thickness of the liquid layer is derived using long-wave theory, from which we also derive a simpler thin-film evolution equation. In the thin-film case, we show that two branches of marginally-yielded static solutions of the evolution equation can be used to both predict the size of initial perturbation required to trigger instability and quantify how increasing the capillary Bingham number (a parameter measuring yield stress relative to surface tension) reduces the final deformation of the layer. Using numerical solutions of the long-wave evolution equation, we quantify how the critical layer thickness required to form a liquid plug in the tube increases as the capillary Bingham number is increased. We discuss the significance of these findings for modelling airway closure in obstructive conditions such as cystic fibrosis, where the mucus layer is often thicker and has a higher yield stress.

Introduction
The surface-tension-driven instability of a liquid layer lining the interior of a cylindrical tube is of physiological importance since, when it occurs in a lung airway, it can cause obstruction or airway closure by redistributing the liquid lining the airway, potentially leading to formation of a liquid plug. The liquid that lines the lung's airways consists primarily of mucus, a non-Newtonian fluid exhibiting various rheological properties such as shear-thinning and viscoelasticity (Hill et al. 2022). Importantly, mucus also has a yield stress, which is significantly increased in diseases such as cystic fibrosis (CF) and chronic obstructive pulmonary disease (COPD) compared to its typical value in healthy lungs (Patarin et al. 2020). Increased prevalence of airway obstruction by mucus plugging is also a key symptom of CF and COPD (Mall 2016). Motivated by this application, we study the effect of viscoplastic liquid rheology on the evolution of a layer coating the interior of a cylindrical tube, which is a simple model for an airway. Additionally, there are numerous potential applications of this class of flow in engineering and industry, as highlighted by Craster & Matar (2009) for thin-film and coating flows, and Balmforth et al. (2014) for thinfilm and free-surface viscoplastic flows.
The surface-tension-driven flow of a viscous film coating a rigid circular cylinder has been well studied in the case that the liquid layer is Newtonian. Goren (1962) identified that a flat layer can be linearly unstable. Everett & Haynes (1972) found and analysed capillarystatic configurations of a volume of liquid inside a tube, which were either annular collars of fluid or liquid plugs. A nonlinear evolution equation using thin-film theory was first derived and solved by Hammond (1983), who found that an initially flat layer evolves into a configuration with large quasi-static annular collars separated by thin films which slowly drain into the collars. Lister et al. (2006) studied the long-time dynamics of the thin-film system, finding that at very long times collars can translate along the tube, potentially consuming other collars in the process, provided other physical effects do not intervene first, while Xu & Jensen (2017) showed how collars can be pinned by wall roughness. Hammond's theory was extended by Gauglitz & Radke (1988) to predict plug formation by retaining certain higher-order terms in the thin-film theory, notably the exact free-surface curvature. This approach provides a composite approximation to the evolution of layers with thickness comparable to the tube radius, which accurately determines capillary-static effects whilst approximating the dynamics well where the layer is thin. They identified a critical average layer thickness, approximately 12% of the tube radius, required for a liquid plug to form during the evolution. A similar composite approximation was compared to full two-dimensional numerical simulations by Johnson et al. (1991); whilst their quasi-one-dimensional theory could predict when a plug would form, it could not capture the genuinely two-dimensional dynamics which occur around coalescence. Otis et al. (1993) derived a similar reduced-order model by making a long-wave assumption when simplifying the governing equations.
At this point, we clarify the distinction between thin-film and long-wave approaches to deriving reduced-order evolution equations: in thin-film theory, it is assumed that the thickness of the layer is much smaller than the radius of the tube, while in long-wave theory, it is assumed that the tube radius is much smaller than the characteristic axial lengthscale of the flow but the layer is not necessarily thin compared to the radius. Making the thin-film assumption results in an evolution equation with the same mobility function as would appear in the planar case, and the curvature of the cylindrical geometry is felt only through the linearised free-surface curvature. In long-wave theory, additional terms appear in the mobility which arise due to the curvature of the geometry, and the full expression for the free-surface curvature is generally retained. Thin-film models benefit from their relative simplicity, but long-wave models capture the effects of the curved geometry more accurately (Camassa & Ogrosky 2015), allowing the dynamics leading to plug formation to be described.
Various physical effects that can modify the evolution of the coating layer have previously been incorporated into models. Halpern et al. (2010) used a long-wave evolution equation to model the effect of viscoelasticity, showing that the critical layer thickness for plug formation is not changed but the time to form a plug can be shortened by increasing the Weissenberg number (a parameter proportional to the relaxation time of the fluid). Romanò et al. (2019) used a volume-of-fluid method to model the pre-and post-coalescence phases of Newtonian plug formation, and recently extended this work to include the effect of viscoelasticity (Romanò et al. 2021). They found that post-coalescence bi-frontal plug growth can induce significant stresses on the tube wall, and that viscoelasticity can induce additional wall stress due to the occurrence of an elastic instability. Erken et al. (2022) studied the instability of a two-layer coating on the interior of a tube as a model for a lung airway which takes into account the periciliary liquid layer that lies beneath the mucus layer and is generally less viscous than mucus. They found that plug formation occurs more quickly in a two-layer model due to the lubricating effect of the base layer, but that the combined critical thickness of the layers required for plug formation can be significantly larger than the single-layer result of Gauglitz & Radke (1988). Halpern & Grotberg (1992) modelled the evolution of a liquid layer coating an elastic tube and subsequently extended the model to include the effect of insoluble surfactants (Halpern & Grotberg 1993). They found that the presence of surfactant can significantly increase the critical layer thickness required to form a plug and can delay plug formation when it does occur, while decreasing the wall stiffness has the opposite effect, decreasing both the critical layer thickness for plug formation and the closure time. Heil et al. (2008) showed that the volume of liquid required for a plug to form is significantly decreased if there is non-axisymmetric collapse of the elastic tube wall. Halpern & Grotberg (2003) developed a thin-film model that included the effect of an oscillating air-flow in the centre of the tube, showing that at certain frequencies of oscillation, air-flow can suppress deformation of the liquid layer. Camassa et al. (2014) developed a long-wave model for gravity-driven flow, and identified families of travelling-wave solutions which they found can be used to predict the critical thickness for plug formation as a function of the Bond number. Camassa et al. (2017) also found travelling-wave solutions for the case of flow driven by air-flow in the centre of the tube and recently Ogrosky (2021) extended the long-wave model to include the combined effects of gravity, air-flow and surfactant.
Turning to viscoplastic flows with applications in airway modelling, Craster & Matar (2000) modelled surfactant-driven flow on a single-layer or two-layer film of viscoplastic or Newtonian fluids, showing that, at least in the single-layer case, yield stress decreases spreading rates and can cause the layer to become frozen in a non-trivial static shape. Modelling of propagation of viscoplastic liquid plugs in tubes and channels has shown that increasing the yield stress increases the stress applied to the wall and increases the thickness of the layer of liquid left behind as a plug propagates (Zamankhan et al. 2012(Zamankhan et al. , 2018. Rupture of viscoplastic liquid plugs has also been modelled both experimentally (Hu et al. 2015) and numerically (Hu et al. 2020), showing that increased yield stress can inhibit plug rupture because a larger pressure drop is required across the plug to make it yield. Recently, Bahrani et al. (2022) proposed a model of elastoviscoplastic plugs, which they validated against experimental results, showing that increased yield stress slows the propagation of a plug but can speed up its rupture since the trailing film thickness is increased. The distribution of mucus throughout a whole lung has also been studied using a viscoplastic model for mucus, showing that the yield stress and the strength of air-flow (in this case modelling air-flow induced by chest physiotherapy) influence the mucus layer thickness in each airway generation (Mauroy et al. 2011(Mauroy et al. , 2015. The exposition of viscoplastic thin-film theory by Balmforth & Craster (1999) has provided the basis for various studies of canonical viscoplastic free-surface flows. In their theory, there are regions of plug-like flow near the free-surface, which have the same structure as the "pseudo-plugs" first identified in a bounded annular flow by Walton & Bittleston (1991). Viscoplastic thin-film theory was used by Balmforth et al. (2000) to study axisymmetrically spreading gravity currents and the work was recently extended to model droplets spreading under surface tension as well as gravity (Jalaal et al. 2021), showing that after spreading, the fluid is frozen in a non-trivial static shape in which the hydrostatic or capillary pressure is balanced by resistance from the yield stress. Jalaal et al. (2021) also compared results for the final shape and size of the droplets computed using thin-film theory to results from computational fluid dynamics (CFD) simulations showing good agreement except when the capillary Bingham number was very large. Gravity-driven flow down inclined planes has been well studied using viscoplastic thin-film theory, as reviewed by Balmforth et al. (2007a). Balmforth et al. (2007b) investigated the surface-tension-driven fingering instability of a film travelling down an inclined plane, finding that increasing the Bingham number (which measures yield stress relative to viscous stress) slows growth of the linear instability and, when it is above a critical value, instability is fully suppressed. Jalaal & Balmforth (2016) also used thin-film theory to model the steady propagation of a bubble through a tube filled with viscoplastic fluid, and compared their results to CFD simulations, showing that thin-film theory predicts the behaviour accurately when the liquid film is thin but less well when the Bingham number is increased and the film is thicker. Viscoplastic flows are often solved using regularisation of the constitutive equation; Frigaard & Nouar (2005) review popular regularisation approaches. Jalaal (2016) introduced a regularisation specifically designed for thin-film flows which we describe in §2.4 and use when solving our evolution equations numerically.
With this study, we aim to quantify the effect of viscoplastic liquid rheology on the surfacetension-driven Rayleigh-Plateau instability of a layer coating the interior of a cylindrical tube. This flow has not previously been studied in the case that the liquid layer is viscoplastic. We ask how the yield stress affects the dynamics during the evolution of the layer and the critical layer thickness required to form a plug. To answer these questions, we derive an evolution equation using long-wave theory, with the detailed flow structure inspired by the viscoplastic thin-film theory of Balmforth & Craster (1999), but we include additional terms arising from the cylindrical geometry which are neglected in the thin-film approximation. We then show how this model reduces, in the appropriate limit, to a thin-film evolution equation, analogous to the Newtonian version derived by Hammond (1983). Other complicating effects are neglected in the model so that the effect of the viscoplastic rheology can be examined in isolation: the tube is rigid, the flow is axisymmetric, surface tension is constant, and the air in the centre of the tube is passive and inviscid. We use the Bingham model for the liquid layer since it is the simplest viscoplastic rheology, without the potentially complicating effects of, for example, shear-thinning, elasticity or thixotropy. We compute marginally-yielded static solutions to the thin-film evolution equation, and show how they can be used to predict both the size of perturbation to a flat layer required to trigger instability, and the final shape of the layer when there is instability. The thin-film theory cannot, however, predict formation of liquid plugs. By solving the long-wave evolution equation numerically, we examine the critical layer thickness required for plug formation and the time taken for plugs to form, quantifying how both can be increased by increasing the capillary Bingham number.
The rest of the paper will be organised as follows. In §2, we formulate the models, presenting the long-wave evolution equation in §2.2, and the thin-film equation in §2.3. A brief discussion of the methods for solving these equations is given in §2.4. Results for thin layers are presented in §3. We discuss a representative numerical solution of the thin-film evolution equation in §3.1, and we examine the behaviour of the layer at long times in §3.2. We compute and analyse static solutions of the thin-film evolution equation in §3.3, and investigate the dependence of the evolution on the initial conditions and the capillary Bingham number in §3.4. Results for layers with finite thickness are presented in §4. We discuss an example numerical solution of the long-wave equation in §4.1 and examine the dependence of the evolution on the capillary Bingham number, layer thickness and initial conditions in §4.2, including discussion of the critical layer thickness for plug formation and the time taken to form a plug. A summary of the results, and a discussion of their significance for modelling lung airways is given in §5.
unyielded yielded yield surface Figure 1: Sketches of the geometry and flow structure in the long-wave model. In non-dimensionalised variables, the free surface is at = and the layer thickness is . The surface = Ψ separates a region of shear-dominated flow near the cylinder wall and a region of plug-like flow near the free surface. We define Ψ ≡ min(1, ) with defined in (2.13) and Ψ = 1 corresponding to regions of unyielded fluid. The thin-film model has qualititavely the same flow structure.
(2.4) The kinematic boundary condition at the free surface is The gas phase applies no shear stress to the liquid film so where are the components of the normal to the free surface, is the constant value of surface tension, and 1 + ( * * ) 2 1 * − * * 1 + ( * * ) 2 (2.7) is the free-surface curvature. At the side boundaries, we impose symmetry boundary conditions, * * = * = * = 0 at = {0, * }.
(2.8) Rather than solve the full Stokes problem defined above, we will derive reduced-order models using long-wave and thin-film theories, presented in §2.2 and §2.3, respectively.
Surface tension at the air-liquid interface introduces an associated energy, proportional to the surface area of the air-liquid interface, * ≡ ∫ * 0 2 * 1 + * * 2 d * . (2.9) In Appendix A.1, we show that the Stokes equations and boundary conditions (2.2)-(2.8) where is the volume of the layer, so the interfacial energy is always decreasing. In the Newtonian problem, the final shape that the layer reaches after its evolution can be found by solving for shapes which locally minimise interfacial energy (Everett & Haynes 1972). However, in the viscoplastic problem, this is not necessarily the case because we expect that the yield stress may freeze some or all of the layer before it has reached a minimal energy state. Analysing the final static shapes of viscoplastic layers will form a large part of our discussion, particularly in the thin-film case (cf. §3.3).

2.2.
The long-wave model We non-dimensionalise the governing equations and boundary conditions (2.2)-(2.8) by defining where hats are used to distinguishˆ ,ˆ ,ˆ andˆ from the scaled, thin-film quantities , , and which we will define in §2.3. After non-dimensionalising, we consider (2.2)-(2.8) in a long-wave limit by introducing a characteristic axial lengthscale for the flow, / , where ≪ 1. There is no assumption at this point that the liquid layer is thin. We then rescale the system using the small aspect ratio, , defininḡ with other variables remaining unstretched. The resulting scaled, dimensionless equations are given in Appendix B.
To derive the long-wave evolution equation, we propose a similar flow structure in the long-wave model to that of the thin-film theory of Balmforth & Craster (1999). Where the fluid is yielded, the flow is separated into a shear-dominated region, 1, adjacent to the no-slip boundary, where the shear stress is large compared with the normal stresses, and a region, < , adjacent to the free surface, where we say the flow is "plug-like" as the axial velocity,ˆ , is independent of (figure 1). We make separate expansions for the velocities and stresses in the shear-dominated and plug-like regions. In the plug-like region, the shear stress is below the yield stress, but the fluid is still yielded because the normal stresses are large enough that the total stress exceeds the yield stress. We solve at leading order in in the shear-dominated and plug-like regions, match the solutions from the two regions together, and finally arrive at the evolution equation which we state below. The full derivation is given in Appendix B.
Following the approach of, e.g., Camassa et al. (2012), we write the evolution equation in terms of the unscaled variables (2.11) rather than the scaled variables (2.12), so does not appear in the equations. However, the limit in which the theory is formally valid remains ≪ 1. From now on, we will use subscripts to denote derivatives. We determine the surface between the shear-dominated and plug-like regions to be is a capillary Bingham number, which measures yield stress relative to capillary stress. We use the hat notation to distinguishˆ from the thin-film version, , which we will introduce in §2.3. The capillary pressure is proportional to the free-surface curvature, and is given bŷ (2.14) We have retained the full expression for capillary pressure (2.14) including all terms which are higher-order in . Although this is not strictly consistent with the asymptotic analysis, it allows the evolution equation to describe capillary static effects accurately and is a widely used device (e.g. Gauglitz & Radke 1988). The axial velocity is defined separately in the shear-dominated and plug-like regions, The function Ψ is defined so that (2.15) applies to the whole layer, including regions of unyielded fluid. Where Ψ = 1, the fluid is unyielded, so there is no motion andˆ = 0. The axial flux,ˆ , is calculated by radially integratingˆ .
The long-wave evolution equation is with non-negative functions 1 and 2 (see figure 9 below) given by The boundary conditions at the sides of the domain are The initial conditions which we impose when solving (2.16) are which corresponds to a flat layer perturbed by a single Fourier mode with wavelength 2 . The constant is the perturbation amplitude and the constant is the ratio of average layer thickness to tube radius when = 0. The constant term in (2.19) is chosen so that the total volume of the layer is independent of for a given .
We derive an expression for the shear stress in Ψ 1 (B 11), which when evaluated at = 1 gives the stress exerted on the tube wall, (2.20) Note that (2.20) only holds in regions where the fluid is yielded (where Ψ < 1). In unyielded regions, the stress is not defined by the constitutive relation (2.3) but we do know that the wall stress must be bounded by the yield stress, so |ˆ | ˆ where Ψ = 1. From (2.9), the dimensionless interfacial energy is (2.21) In Appendix A.2, we deduce directly from (2.14)-(2.18) that ˆ 0. Hence, the result (2.10) is preserved in the long-wave theory.

The thin-film model
We now derive the analogous evolution equation for a thin film, no longer requiring ≪ 1, but assuming |1− | ≪ 1. The thin-film evolution equation can be derived using the approach of Balmforth & Craster (1999), but here we derive it by taking a thin-film limit of the longwave system (2.13)-(2.18). The thin-film approximation acts to flatten the geometry, so that terms in the long-wave mobility (2.17) which arise due to the curvature of the geometry are negligible. The effect of the cylindrical geometry is then only felt through the free-surface curvature, which is linearised.
We consider a layer with characteristic thickness , and now let ≪ 1. We rescale time, defining ≡ 3ˆ , then define the dimensionless film thickness, ( , ), which satisfies It is convenient to define a radial coordinate, , measured from the no-slip boundary, which satisfies = 1 − . Then the free surface is located at = . Substituting (2.22) into (2.14) givesˆ = −ˆ = −1 − ( + ) + ( 2 ). We define the thin-film curvature and capillary pressure as ≡ (ˆ − 1)/ and ≡ (1 +ˆ )/ , then linearise in , so that the pressure gradient driving the flow is We define the thin-film capillary Bingham number as As before, we augment this definition so that it holds in yielded and unyielded regions: we define ≡ max(0, Y), which obeys Ψ = 1 − + ( 2 ), with = 0 corresponding to regions of unyielded fluid. Where > 0, the fluid is yielded, and the flow structure is qualitatively the same as in the long-wave model. The value of then indicates the boundary between the shear-dominated and plug-like regions of flow. The axial velocity (2.15) becomeŝ = 3 + ( 4 ), where the thin-film axial velocity is . (2.26) Finally, substituting (2.22)-(2.25) into (2.16)-(2.17), and linearising in , gives Equation (2.27), with definitions (2.23) and (2.25), is the thin-film evolution equation. In the thin-film limit, the boundary conditions (2.18) become where we define the thin-film flux as ≡ 2 ( − 3 )/6. Note that in the Newtonian problem (e.g. Hammond 1983), enforcing zero flux at = {0, } is equivalent to enforcing zero third derivative, = 0 or = 0. Here, the zero flux conditions (2.18) and (2.28) are preferable because the third derivatives, and , generally become discontinuous at = {0, } during the evolution. This also occurs at any interior points where the direction of flow changes (cf. figure 2c). This is an inconsistency in the theory which could be resolved by finding a solution in the inner region (likely of axial length ( ) in the thin-film system or ( ) in the long-wave system) around each of these points and matching these to the global outer solution which we compute. Following the approach of, e.g., Balmforth et al. (2000), we do not solve in these inner regions and assume that the solution which we compute captures the global dynamics of the layer sufficiently accurately.
After linearising in and combining with (2.22), the initial condition (2.19) becomes which is the initial condition we will use when solving (2.27). We define the thin-film wall shear stress, ≡ˆ / 2 , then linearise in to get where we used (2.25) in the second equality. Note that (2.30) only holds in regions where the fluid is yielded ( > 0), but we have the bound | | in unyielded regions ( = 0). The interfacial energy (2.21), when expanded in powers of , becomes where 0 is the (constant) total volume of the layer. We show in Appendix A.3 that (2.23), (2.27) and (2.28) imply 0 for ≪ 1. Hence, the thin-film approximation also preserves the result (2.10).

Solution methods
When solving both the thin-film and long-wave equations, we choose the domain length to be = √ 2 . This length corresponds to the half-wavelength of the most unstable mode in the Newtonian linear stability analysis (Hammond 1983). The instability in the viscoplastic problem is inherently nonlinear, but this choice for allows direct comparison to previous literature on the Newtonian and viscoelastic versions of the problem (e.g. Gauglitz & Radke 1988;Halpern et al. 2010). We found that small changes in do not qualitatively affect our results so = √ 2 can be considered a representative domain length. The form of perturbation in the initial conditions, (2.19) or (2.29), then corresponds to the single unstable Fourier mode that exists in the domain. In the long-wave theory, is defined as the ratio of tube radius to a typical axial lengthscale. If that axial lenthscale is taken to be the wavelength of the initial disturbance, then = 1/(2 ) = 1/(2 √ 2 ). Shorter wavelength structures also develop in the thin-film and long-wave simulations (cf. capillary waves discussed in §3.1), testing the validity of the long-wave theory. Pending validation by computations of the full problem (2.1)-(2.8), we anticipate that our results provide a good approximation to the true behaviour, with additional accuracy gained from retaining the exact expression forˆ in (2.14). We do not observe any significantly different behaviour when applying periodic boundary conditions at the sides of the domain compared to the boundary conditions (2.18) or (2.28), validating our use of the latter in all the results presented.
When solving the systems (2.13)-(2.18) or (2.23)-(2.28) numerically, we use the method proposed by Jalaal (2016). The evolution equations are regularised by redefining = 10 −6 after confirming this is small enough that the results are not sensitive to the precise value of . For example, for the simulation presented in figure 6 below, the absolute errors in max ( , = 120) and (the time to form a plug) are bounded above by 800 2 and 40000 2 , respectively, for all 10 −6 10 −3 , which we find to be typical of convergence rates in simulations. Where = or Ψ = Ψ , there is a very weak regularisation-induced-flow, but since is chosen small enough for this flow to be negligible, we treat these regions as unyielded, treating = or Ψ = Ψ as equivalent to = 0 or Ψ = 1. The regularised equations are solved using the method of lines: the spatial derivatives are approximated using second-order centred finite differences and the resulting system of ODEs is solved through time using a stiff solver in M . We have confirmed that the number of spatial grid points used is large enough that the precise value does not affect our results.

Results: thin-film theory
3.1. Time evolution of a thin layer Figure 2(a) shows snapshots from a numerical solution of the thin-film equations (2.23)-(2.29) with = 0.05 and = 0.2. For any viscoplastic simulation, the initial perturbation must be sufficiently large in order to trigger instability because a sufficiently large pressure gradient must be created to overcome the yield stress and make the fluid yield. Here, is large enough that the fluid in the centre of the domain yields, but = 0 near the boundaries so the fluid there is initially rigid (figure 2a, = 0). There is an initial period during which there is minimal deformation in , but deforms significantly and, by = 20, > 0 for all ∈ (0, ) indicating that the whole layer is yielded. This initial yielding period causes a delay in the growth of the instability, as can be seen when max ( , ) is compared to the same quantity from a Newtonian ( = 0) simulation in figure 2(b). After the initial yielding period, there is a period of significant deformation of the free surface. This coincides with a peak in max ( , ) and a peak in the wall shear stress, | |, at around = 100. This indicates that there is significant shear during this period, with more of the layer exhibiting shear-dominated flow and the region of plug-like flow becoming smaller.
At late times, the layer relaxes slowly towards a final marginally-yielded static shape in which → 0 across the whole layer (figure 2(a), = 600). Figure 2(b) shows that max decays towards zero at a rate proportional to −1 . From (2.30), we note that as → 0, | | → across the entire layer, which indicates that even when the layer reaches its final static shape, viscous and capillary effects apply a uniform stress on the tube wall equal to the yield stress. Figure 2(b) also indicates that max | | decays towards at a rate proportional to −1 , in contrast to the Newtonian result that the peak wall shear stress decays towards zero at a rate proportional to −1/4 (Jones & Wilson 1978; Hammond 1983). The late-time shape of the layer consists of a large collar of fluid around = and a small collar around = 0. Figure 2(b) shows that the final value of max is lower for the = 0.05 solution compared to the Newtonian solution. This is because, unlike in the Newtonian evolution, not all of the fluid drains into the large collar at late times. Instead, some is trapped in the small collar so the peak height of the layer is decreased. Thus, the yield stress inhibits the growth of the instability. In §3.2 and §3.3, we quantify how increasing affects the final marginally-yielded static shape of the layer.
In the early-time period of gradual yielding, capillary-wave-like structures form ahead of the travelling yield surfaces (where makes contact with zero). They can be identified by observing the structure of and the pressure gradient ( figure 2c). There is a jump discontinuity and a change of sign in between each of these structures, indicating that the direction of flow reverses. At each point that passes through zero, we also have = 0. We expect there to exist additional waves with smaller wavelengths ahead of those observed in figure 2(c), but the numerical method can only resolve the largest few since it is limited by the size of the finite difference grid spacing. The observed structures resemble closely the capillary waves identified by Balmforth et al. (2007b) and Jalaal & Balmforth (2016), and studied in detail by Jalaal et al. (2021) in the context of spreading viscoplastic droplets. They are a feature common to surface-tension-driven viscoplastic flows, occurring when a yield surface advances into a region of unyielded fluid. Unlike in previously studied flows, here the capillary waves, in general, only exist transiently, during the early time period of gradual yielding until the whole layer yields. However, we find that for some values of and (mostly very large ), one or more of the jump discontinuities in which develop can persist for the whole evolution. In §3.4, we discuss how this phenomenon can affect the final static shape of the layer, and present criteria for it to occur. Until then, we focus on the case that the flow is unidirectional ( > 0 for all ∈ (0, )) after the early time capillary waves have passed.

Late time asymptotics for the thin-film evolution equation
To analyse the late time dynamics of the layer, we look for a solution in which = Y → 0 as → ∞. Numerical simulations suggest that decays like −1 (figure 2b), so we make the expansions The −1 rate of decay of towards a steady state is also consistent with the numerical results in figure 2. For this analysis, we assume the capillary pressure is monotonic, or equivalently the pressure gradient is one-signed, − 0, − 0, < 0 for all ∈ [0, ]. This is equivalent to assuming unidirectional flow at late times. Substituting the expansions (3.1) into the definition of Y (2.25) and the evolution equation (2.27) gives The boundary conditions (2.28) imply is a static solution of the evolution equation (2.27) in which Y = 0 uniformly. This is not a capillary-static solution in which the pressure is everywhere uniform; instead it is a state in which the layer is uniformly marginally yielded. From (2.30) we note that, in the static solution, = uniformly, indicating that there is a stress being applied in the positive -direction, but it is resisted by the yield stress, preventing flow. The functions 1 ( ; ) and 1 ( ; ) quantify the rate at which the layer approaches the static solution at late times.
The equations (3.2a)-(3.2c) were solved subject to (3.3) and (3.4) using a boundary value problem solver in M . A solution for = 0.05 is shown in figure 3 with comparison to the final snapshot of the numerical simulation from §3.1. Figure 3(a) shows the agreement is very good between 0 and the late-time shape of the layer from the numerical solution. Figures 3(b) and 3(c) show that 1 and 1 approximate well the rates of decay of towards 0 , and towards zero, respectively. This indicates the expansion (3.1) accurately describes the late-time dynamics of the evolution and confirms the ( −1 ) rate of decay in determined empirically in Fig 2(b).

Static solutions
The marginally-yielded static solutions, 0 ( ; ), can predict the final state of the layer. To investigate the dependence of 0 on , we solve (3.2a) with (3.3) and (3.4a), varying . We find that there exists a value * ≈ 0.163 such that for all 0 < * exactly two solutions for 0 exist, and for > * no solutions exist. There is a bifurcation at = * . Figure  4(a) shows max 0 ( ; ) for all the solutions, which always coincides with 0 ( = ; ). to max 0 ( ; * ) ≈ 1.98. This indicates that increased yield stress can significantly inhibit deformation of the film.
The lower-branch solutions are near-flat for small , becoming more deformed as is increased. In numerical simulations, an unstable layer evolves away from a near-flat configuration towards a strongly deformed (upper-branch) static shape. This suggests that the upper-branch solutions are stable and the lower-branch solutions are unstable. We confirm the stability of the two branches using asymptotic analysis near to the bifurcation, ≈ * , presented in full in Appendix C. We find that ( , ) ∼ 0 ( ; * ) + A ( ) 1 ( ) as ≡ √ * − → 0, where 1 ( ) is a solution to a linear ODE, = 3 is a slow timescale, and A ( ) is an amplitude function which solves an ODE of the form where 0 and A 0 are positive constants. Equation (3.5) has two fixed points, A = ±A 0 , and we compute A 0 ≈ 2.20. Solutions evolve away from the negative fixed point towards the positive one (figure 10b below), indicating that the positive fixed point is stable and the negative one unstable. These fixed points correspond to two static solutions for which are the upper-and lower-branch solutions, respectively. With this stability result, we identify the bifurcation at = * as a saddle-node bifurcation. We also approximate the location of the branches in figure 4(a) by Since the lower-branch solutions satisfy Y = 0, they are marginal states between rigid layers (Y 0) and fully yielded layers (Y > 0). We expect that if a layer is initially more deformed than the lower-branch solution it will be yielded and unstable, but if it is less deformed initially it is likely to be rigid and thus stabilised. We provide evidence from numerical simulations to confirm this in §3.4, where we show that the lower-branch static solutions correspond almost exactly to the minimum amplitude of initial perturbation required to trigger unstable growth.
Asymptotic analysis for small shows that the lower-branch solutions have the regular expansion, Taking the maximum value of (3.7) gives which approximates the lower branch in figure 4(a). To approximate the upper-branch solutions for small , we construct a solution by matched asymptotic expansions. The analysis is presented in full in Appendix D and illustrated in figure 11 below. In addition to the large collar around = and the small collar around = 0, we identify a third, inner region located around = − where 0 ∼ ( 2 ). In contrast to the Newtonian problem, where the inner region is described by an ODE of the form 3 = 1 (Hammond 1983), here the relevant ODE is of the form = 1 (see D 5). The difference arises because in the Newtonian problem there is constant flux across the inner region during the late-time draining regime, while here there must be constant stress across the inner region since 0 is marginally yielded. After expanding and solving for 0 ( , ) in each of the three regions and matching the solutions, a composite approximation to 0 ( , ) is found. This also provides an approximation to the maximum value, where 1 is a constant which depends on . For = √ 2 , we compute 1 ≈ −0.105. Equation (3.9) approximates the upper branch in figure 4(a).
At the saddle-node bifurcation, = * , the two static solutions annihilate each other, so it is only possible for the layer to select the upper-branch solution if < * . Note that in this section, we have only computed the static solutions that have monotonic pressure, but other static solutions may exist. We show in the next section that even when < * , the layer may select a different final static shape, depending on the initial conditions of the layer, and that it is possible to have some yielding and unstable growth for some > * if the initial perturbation is sufficiently large.

Dependence on capillary Bingham number and initial conditions
The final shape of the layer can depend on its initial conditions as well as . To investigate this dependence, we solve the thin-film equations (2.23)-(2.29) numerically for a range of values of and a range of initial perturbation amplitudes . We run simulations on a regularly spaced grid of points in the range (0 0.5, 0 0.99). All simulations are run to a fixed, long time, which we choose to be = 1000. In figure 5(a), the final maximum height, max ( , = 1000), is plotted for each simulation. Figure 5(a) shows that there are three qualitatively different possible outcomes for an evolving layer, depending on the values of and . First, there is a region for small and large with max ( , = 1000) = 1 + , so the final maximum height is equal to the initial maximum height. In these cases, the initial perturbation does not generate a large enough pressure gradient to make the layer yield, so the yield stress entirely suppresses unstable growth. This region also extends to all of > 0.5 for all 0 < 1. We call this the yield-stabilised region. We can seek a bound for the yield-stabilised region using the minimum amplitude, = ( ), such that the whole layer is initially unyielded if and only if < . We identify that is the value of such that the initial condition (2.29) makes  1000). The data are linearly interpolated to produce the black contour lines, which are evenly spaced. When is small and is large, the layer is yield-stabilised: no unstable growth occurs. The critical amplitude for any yielding to occur, ( ) (dashed red), defined in (3.10), is a strict lower bound on the boundary of the yield-stabilsed region. When there is growth, the final shape either has monotonic or non-monotonic pressure, depending on and . The quantity 1 − 0 ( = 0, ) (magenta) from the static solutions (figure 4) predicts the boundary of the monotonic pressure region. The maximum value of along the magenta curve is * . Two large black dots indicate the location of the example solutions shown in (b) and (c), with , 70 (scaled for clarity) and plotted at = 1000. (b) shows a solution in the monotonic pressure region ( 0) and (c) shows a solution in the non-monotonic pressure region ( changes sign once in the domain).
Y non-negative at exactly one point in the domain, and thus we find is given implicitly by where = / = √ 2/2. The curve = ( ) provides a strict lower bound on the boundary of the yield-stabilised region (figure 5a). A small number of yield-stabilised simulations have > ( ): in these simulations the fluid initially yields a small amount in the centre of the domain but rigidifies before the sides of the domain yield, so there is no growth in max .
In the second possible outcome, there is unstable growth and the layer evolves towards the upper-branch static solution with monotonic pressure which we computed in §3.3. This occurs for a large set of and as indicated in figure 5(a), with max ( , = 1000) being independent of in this region since all simulations approach the same final shape for a given . Figure 5(b) shows the final shape for a simulation with = 0.4, = 0.1, which is within the monotonic pressure region. The pressure gradient is non-positive meaning the flow is unidirectional in the positive -direction at = 1000 when the simulation is stopped.
is positive but close to zero for all ∈ (0, ) at = 1000, so the layer is almost rigid and is very close to the upper-branch static solution. Within the monotonic pressure region in figure 5(a), the value of max ( , = 1000) decreases as is increased, consistent with the decrease shown in the upper branch in figure 4(a).
The final possible outcome involves unstable growth of the layer leading to the final static shape having non-monotonic pressure. Figure 5(a) shows that this occurs mainly when is very large. In these simulations, max ( , = 1000) depends on both and and the final shape is not predicted by the upper-branch static solutions in figure 4(a). An example of a late-time shape with non-monotonic pressure is shown in figure 5(c) from a simulation with = 0.95, = 0.1. Comparing this with figure 5(b), there is less fluid in the collar near = 0 and so the collar near = is larger, even though is the same. In figure 5(c) there is exactly one point in the domain where the sign of changes, which corresponds to a point where the direction of flow changes. For most values of and in the non-monotonic pressure region, the final shape selected by the layer has exactly one sign-change in but for very close to 1 we found that there can be more. The sign-changes in develop during the early-time yielding period, caused by the presence of capillary waves near the yield surfaces (figure 2c). Most of the sign-changes exist only during this early-time period, but we see here that one or more can persist and affect the late-time dynamics, suggesting that in these cases the early-time capillary waves can influence the entire evolution. In general, the location of the sign-change(s) in affects the final shape of , and the location of the sign-change(s) depends on the initial shape of the layer.
There is a sharp boundary in the numerical data between the yield-stabilised and monotonic pressure regions, across which max ( , = 1000) jumps significantly. There is also a clear boundary between the monotonic and non-monotonic pressure regions, indicated by where the contours of max ( , = 1000) begin to curve. Figure 5(a) shows that we can predict the locations of both boundaries using the quantity 1 − 0 ( = 0; ), where 0 ( ; ) are the static solutions with monotonic pressure computed in §3.3.
First, we give an explanation for why 1 − 0 ( = 0; ) coincides with the boundary between the yield-stabilised and monotonic pressure regions. Consider a simulation with just above the boundary, so the fluid yields just enough to trigger instability (e.g. figure 2). After the initial period of gradual yielding, Y is very close to zero but positive everywhere. At this point the layer's shape is very close to the lower-branch static solution with the same , which has Y = 0 everywhere. Since the fluid at = 0 is unyielded for most of the initial period of the evolution in which the layer gradually yields, the height at = 0 does not change significantly in this period. Hence, the height at = 0 must have initially been very close to 0 ( = 0; ), the height of the lower-branch static solution at that point. The initial height at = 0 is (0, 0) = 1 − , so the boundary between the yield-stabilised and monotonic pressure regions is predicted by = 1 − 0 ( = 0; ), and the lower-branch static solutions effectively correspond to the minimum amplitude of perturbation required to trigger instability.
Secondly, we give an explanation for why 1 − 0 ( = 0; ) coincides with the boundary between the monotonic and non-monotonic pressure regions in figure 5(a). To do this, we determine a condition for a final solution with monotonic or non-monotonic pressure to be selected during the evolution. If the initial height of the layer at = 0 is smaller than the height of the corresponding (i.e. for the same ) upper-branch static solution at = 0, then there must be fluid flow towards = 0 in the negative -direction during the evolution if this solution is to be selected. This would mean that the flow at late times must not be unidirectional, and hence the pressure of the final shape must be non-monotonic. (The transient early time capillary waves in these simulations create flow reversal but it is very weak so can be neglected in this argument.) So, if ( = 0, = 0) < 0 ( = 0; ), for a The simulation is stopped when max = 0.7/ as it is clear that a plug will form. The evolution of max and max |˜ | for a Newtonian ( = 0) simulation with the same and shows plug formation occurring earlier, ≈ 70. (c) Snapshots of the = 0.05 simulation at evenly spaced time-points between = 50 (darkest lines) and = 130 (lightest lines). Inset showsˆ becoming equal to zero across progressively more of the domain, indicating that the small collar of fluid around = 0 is rigidifying. There is a small jump in whereˆ makes contact with zero, but stillˆ 0 everywhere.
given value of , then the layer will select a final shape with non-monotonic pressure. Noting again that (0, 0) = 1 − , we see that > 1 − 0 ( = 0; ) is an equivalent condition for the layer to select a final shape with non-monotonic pressure. The results in figure 5(a) are specific to the sinusoidal form of initial perturbation used (2.29). However, we also ran simulations with initial conditions of the form ( , 0) = 1 + tanh (2 − ) and the results for max ( , = 1000) were qualitatively, and largely quantitatively, the same. The quantity 1 − ( = 0; ) was still found to accurately bound the monotonic pressure region.
In this section, we have illustrated the complex dependence of the evolution of a thin layer on and , and shown how the static solutions with monotonic curvature can provide insight into the dynamics of the layer. However, thin-film theory cannot capture the full range of possible dynamics for the system because the volume of fluid in a thin layer is too small to form a liquid plug in the tube. We now address this by using long-wave theory to model layers with finite thickness.

Time evolution of a layer with finite thickness
Figure 6(a) shows snapshots from a numerical solution of the long-wave equations (2.13)-(2.19) with film thickness = 0.125, capillary Bingham number = 0.05 and initial perturbation amplitude = 0.2. For ease of comparison with the thin-film results, we describe solutions in terms of thin-film parameters and variables: instead ofˆ we use =ˆ / 2 , instead ofˆ we use = 3ˆ , instead of ( ,ˆ ) we use ( , ), instead of Ψ( ,ˆ ) we useˆ ( , ) ≡ (1 − Ψ)/ , and instead ofˆ we use˜ ≡ˆ / 2 . The early-time behaviour is qualitatively the same as in the thin-film simulations. There is a delay to the growth as the fluid gradually yields (figure 6b) and capillary waves develop, which are qualitatively the same as in the thin-film case. There is then a peak in max ˆ at around ≈ 30, coinciding with significant deformation of the layer. Figure 6(b) also shows that there is a small associated peak in the wall shear stress, max |˜ |, which occurs slightly later, around ≈ 40. The layer then evolves towards a shape with a large collar near = and a smaller collar near = 0 (figure 6a, = 70) Gauglitz & Radke (1988) identified the critical thickness required to form a plug in their Newtonian simulations as = 0.12. Since = 0.125 is larger than this critical value, we expect plug formation may be possible in this simulation. Indeed, figure 6(b) shows that at around = 140, max begins to rapidly increase towards the centre of the tube, which is located at 1/ = 8. The long-wave theory cannot model coalescence so, following the approach of, e.g., Halpern et al. (2010), we stop simulations when max = 0.7/ , but when we do run the simulation further max rapidly approaches 1/ , so it is clear that a plug will form. We denote the time taken to form a plug as , and use the time at which we stop the simulation as a proxy for . Figure 6(b) shows that the plugging time, ≈ 140, is significantly longer than the plugging time for a Newtonian ( = 0) simulation, ≈ 70. This is partly due to the delay caused by the initial yielding period, and partly because the rheology slows down the subsequent period of growth. Throughout the evolution, max |˜ | is larger for the = 0 simulation than for = 0.05, suggesting that the yield stress decreases the wall shear stress during this pre-coalescence phase of plug formation. Note that max ˆ and max |˜ | increase rapidly around ≈ , suggesting that the fluid in the large lobe is strongly yielded and the wall shear stress increases as plug formation occurs. However, we cannot expect the theory to remain accurate during this period since the assumption that radial velocities are weak no longer holds. As in the Newtonian problem (Johnson et al. 1991), fully two-dimensional theory is required to capture the coalescence phase of the evolution.
A new phenomenon which we have observed only in long-wave simulations is that the small collar of fluid which forms near = 0 can rigidify during the evolution. At = 30 (figure 6a), the layer is fully yielded withˆ > 0 for all ∈ (0, ). Figure 6(c) shows that a yield surface (whereˆ makes contact with zero) then travels from = 0 at around = 50 to around = 1.8 at = 130, indicating that almost the entire small collar has rigidified by this point. The speed at which the yield surface travels through the domain decreases slightly through the evolution. We observe a small jump inˆ at the the point whereˆ becomes zero, butˆ remains non-positive everywhere including in the rigid region. We generally observe this rigidification of the small collar in our simulations whenever there is enough fluid to form a plug, i.e. 0.12. For thinner layers, 0.11, we generally do not observe this phenomenon; instead,ˆ decays to zero from above everywhere in the domain, qualitatively the same behaviour as we observed in thin-film simulations (e.g. figure 2).

Dependence on capillary Bingham number, layer thickness and initial conditions
We investigate the dependence of the long-wave evolution on the capillary Bingham number, layer thickness, and initial conditions by running large numbers of numerical simulations varying , and . This allows us to examine the dependence of the critical thickness required to form a plug, , on and . In figures 7 and 8, we plot the final maximum heights of the layers from the numerical simulations. We run all simulations to = 1000, or if a plug begins to form before this time, the simulation is stopped and the stopping time is identified as .
In figure 7, we vary and between simulations while the initial perturbation amplitude is fixed at = 0.25. There are three distinct regions in the data, corresponding to three qualitatively different outcomes for the layer. Firstly, when is sufficiently large, there is no growth so we say the layer is yield-stabilised. As in the thin-film case in §3.4, we can find the minimum amplitude for yielding, = ( , ), such that the layer is initially fully rigid if and only if . This defines a corresponding capillary Bingham number, = ( , ), such that for a given and , the layer is fully rigid if and only if . We compute ( , ) numerically by finding = such that the initial condition (2.19) makes = 1 at exactly one point in the domain, where is defined in (2.13). Figure 7 shows that = ( , = 0.25) provides a strict upper bound on the boundary of the yield-stabilised region. The second possible outcome for the layer is that there is unstable growth but a liquid plug does not form, and instead the final shape at = 1000 is a two-collar configuration, like in the thin-film simulations. Figure 7 shows that this occurs for roughly 0 0.1 and 0.12, and the figure also indicates how the final peak height of the large collar, max ( , = 1000), depends on both and .
In the third possible outcome, the layer forms a liquid plug, which occurs when the layer is sufficiently thick, and is not too large for it to be yield-stabilised. The boundary of the plugging region in figure 7 corresponds to the critical thickness required for plug formation to occur, , which can be seen to depend strongly on . For very small , the critical value is around ≈ 0.12, which is consistent with the results of Gauglitz & Radke (1988). As is increased, first increases slowly, to somewhere in the range 0.125 < < 0.1275 when = 0.11, with two-collar final shapes being observed at = 0.125 when = 0.11. ( , = 1000). The data are interpolated linearly to produce the black contour lines, which are evenly spaced. Grey points correspond to simulations which are stopped early due to a plug forming. In (c), grey contours in the plugging region indicate the plugging time, . The critical amplitude for any yielding to occur, ( , ) (dashed red), provides a strict lower bound on the boundary of the yield-stabilised region. Unlike in the thin-film case (figure 5a), the quantity 1 − 0 ( = 0; , ) (solid magenta), where 0 are the static solutions computed in Appendix E, does not provide useful information on the final shape of these layers. This suggests that for a few values of around ≈ 0.125, the layer can exhibit plugging, two-collar or yield-stabilised behaviour, depending on the value of . For 0.12, the plugging region is bounded by the yield-stabilised region and increases rapidly as is increased. Figure 7 also shows how the plugging time, , decreases as is increased and increases as is increased. There is a rapid increase in near to the boundary of the plugging region, suggesting that we have located the boundary accurately by running simulations to = 1000. Simulations in the plugging region which are near the boundary spend a long time in a near-static two-collar shape (e.g. figure 6a, = 70) before eventually transitioning to form a plug.
In figure 8 we investigate the dependence of the evolution on and , for = 0.12, 0.125, 0.13. We can again identify a yield-stabilised region, a two-collar region, and a plugging region for each value of . The boundary of the yield-stabilised region does not change significantly between = 0.12 and = 0.13, which is consistent with the results in figure 7. This boundary corresponds to the minimum amplitude of perturbation required to trigger instability, which can be seen to strongly depend on . Again, we put a strict lower bound on the boundary of this region using ( , ), the minimum amplitude for any yielding to occur. The size of the two-collar region decreases quickly as is increased, and it has almost entirely disappeared when = 0.13. This suggests that for 0.13, as long as is large enough to trigger growth, a plug is guaranteed to form. When the two-collar region does exist ( figure 8a,b), the location of its boundary with the plugging region depends on and . When is small, this boundary is largely independent of , but when is large, the plugging region extends to higher . We propose that this is because highly deformed initial conditions place most of the fluid near = , so, compared to simulations with small , less fluid is trapped in the small collar near = 0, making the large collar larger and more unstable to plug formation. The results in figure 8 show that the boundary of the plugging region depends on as well as on and . Hence, the critical thickness for plug formation, , must also depend on as well as . The plugging time also depends on both and , as indicated in figure 8(c), with increasing when is increased and decreasing when is increased. As in figure 7, increases rapidly near to the boundary of the plugging region.
In the thin-film problem, the quantity 1 − 0 ( = 0; ) from the static solutions could be used to predict the final shape of the layer in a large number of cases (figure 5a). We have also computed static solutions, 0 ( ; , ), for the long-wave problem by solving = 1 and assuming monotonic pressure (see Appendix E). The quantity 1 − 0 ( = 0; , ) is plotted for = 0.12, 0.125, 0.13 in figure 8. It predicts the threshold amplitude required to trigger instability for small , but underestimates it for larger , and the prediction becomes less accurate as is increased. The upper branch of the curve does not appear to be correlated with the numerical results in any way. In numerical simulations with 0.12, we generally observe that rigidification of the small collar near = 0 (figure 6b) occurs whether a plug forms or not. When this rigidification happens, the final static shape of the layer does not satisfy = 1 everywhere, so the layer selects different static solutions than the ones we have computed. The static solutions we have computed may predict the final shapes of the layer when is smaller, but they cannot do so when there is rigidification of the small collar or when there is plug formation.
The results in figures 7 and 8 are specific to the sinusoidal form of initial perturbation used (2.19). However, we also ran simulations with initial conditions of the form ( , 0) = 0 + tanh (2 − ), where 0 is a constant chosen so that the total volume of fluid is independent of , and found the same qualitative behaviour and only minor quantitative differences. This suggests that the observed behaviour is not strongly dependent on the exact form of initial perturbation.

Discussion
To summarise, we have quantified how viscoplastic rheology can either inhibit growth of, or fully suppress, the surface-tension-driven instability of a layer of liquid coating the interior of a cylindrical tube. We found that for both thin layers and layers with finite thickness, the final shape after evolution depends sensitively on the capillary Bingham number, , as well as on the initial amplitude of perturbation, . Using thin-film theory, we showed that when is below a critical value, which depends on , there is no unstable growth because the fluid does not yield. When there is unstable growth, the final shape of the layer either coincides with the marginally-yielded static solution, 0 ( ; ), from the upper branch of figure 4(a), or the final shape has non-monotonic pressure. Figure 5(a) shows that the quantity 1 − 0 ( = 0; ) from the static solutions accurately predicts both the minimum required to trigger instability, and the large set of and for which the final shape of the layer is 0 ( ; ). By solving the long-wave evolution, we quantified how the critical layer thickness, , required to form a liquid plug is increased by increasing . Figure 7 shows that can be increased significantly beyond the Newtonian value of ≈ 0.12 found by Gauglitz & Radke (1988), primarily because, when is sufficiently large, there is no yielding so no unstable growth. For 0.12 0.13, it is possible for there to be no unstable growth, unstable growth leading to plug formation, or unstable growth with no plug formation, depending on the values of and (figure 8). For > 0.13, if is large enough to trigger instability, a plug will form.
One application of our results is to modelling mucus flow and airway closure in lungs. We used a thin-film capillary Bingham number, ≡ / 2 , to measure the relative strength of the yield stress in the flow. To estimate for a 12 th generation healthy airway, we propose the following typical values: airway radius = 0.4mm (Hsia et al. 2016), surface tension = 30mNm (Chen et al. 2019), and mucus yield stress = 0.27Pa (Patarin et al. 2020). For a mucus layer with thickness = 0.125, this gives ≈ 0.2. Figure 8(b) shows that for = 0.2 we expect plug formation to be possible, but only for 0.6, i.e. only if the mucus layer is significantly deformed initially. Patarin et al. (2020) measured the yield stress in cystic fibrosis (CF) mucus to be 6.34Pa, which would correspond to ≈ 5 when = 0.125. Figure 8(b) shows that = 5 is well inside the yield-stabilised region for all , suggesting that airway closure would not occur via this mechanism for these parameter values. However, other key symptoms of CF are increased volume of mucus in airways and surfactant deficiency (Tiddens et al. 2010), which would correspond to increases in and and, hence, a potentially significant decrease in , making plug formation more likely to occur. Thus, the net effect of CF symptoms on the likelihood of airway closure by this mechanism is not obvious. Experiments and numerical modelling of plug rupture (Hu et al. 2015(Hu et al. , 2020 suggest that, once an airway does close, increased yield stress makes airway reopening more difficult, which would contribute to the increased prevalence of plugged airways in CF. Our results also suggest that airway closure could be triggered if is suddenly decreased, which could be caused by applying certain therapies which are commonly used in CF, such as mucolytics which decrease the mucus yield stress (Patarin et al. 2020) or expectorants which increase the volume of liquid (Donaldson et al. 2006). However, detailed modelling of the effect of such therapies would be required to confirm this conjecture.
We have also shown that yield stress can delay plug formation when it does occur (figures 6c, 7 and 8c). If we set = 10 −2 Pa s then the dimensional plugging time is * = ( / 3 ) s ≈ 0.07 s for = 0.125. For the Newtonian simulation in figure 6(c), ≈ 70, corresponding to * ≈ 5 s, which is approximately the length of a breathing cycle. The = 0.05 simulation in figure 6(c) takes about twice as long to form a plug, so * is likely to be longer than a breathing cycle meaning airway closure is less likely to occur via this mechanism. If we relate to the measured viscosity of mucus (Lai et al. 2009), = 10 −2 Pa s is a feasible value but it could also be significantly larger, meaning plug formation for these simulations could be on the scale of minutes or hours instead of seconds. The layer thickness also strongly influences (figure 7), and also * depends inversely on 3 , so a modest increase in layer thickness can significantly decrease the time taken for plug formation to occur.
Our results also suggest that the shear stress exerted on the tube wall during the precoalescence phase of plug formation can be decreased by yield stress (figure 6b). This has physiological significance because a large shear stress exerted on an airway wall may cause epithelial cell damage (Huh et al. 2007). However, we expect that the wall shear stress is likely to be much larger in the post-coalescence phase, as is the case when the liquid is Newtonian (Romanò et al. 2019), so we cannot make conclusions about the effect of yield stress on wall shear stress during the entire closure process. Additionally, when the layer is too thin to form a plug, the wall shear stress at late times approaches the yield stress (figure 2b), so in these cases it increases as yield stress is increased.
We have focused on investigating the effect of viscoplastic rheology, so other physical effects which are relevant to airway modelling have been neglected. Various extensions to our work could be made to investigate how viscoplastic effects interact with, for example, shear stress induced by air-flow, elastic tube walls or surfactant, all of which have been studied in the case that the liquid is Newtonian (Halpern & Grotberg 1992, 1993, 2003. Additionally, in order to isolate the effects of the viscoplastic rheology, we have not incorporated shearthinning or viscoelastic rheologies, which are also known to be exhibited by mucus (Hill et al. 2022). It remains an interesting open question how the addition of other rheological properties would affect the dynamics of a viscoplastic layer as studied here.
There are some limitations to the thin-film and long-wave theories that we have used to derive reduced-order models. Thin-film theory cannot predict the formation of liquid plugs, and the quasi-one-dimensional long-wave theory cannot capture the fully two-dimensional dynamics which develop as a liquid plug is forming (requiring simulations to be stopped just before coalescence). Viscoplastic thin-film theory is known to break down at points where the direction of flow changes and the pressure gradient has a jump discontinuity (Balmforth et al. 2000), and we observed this same behaviour in our long-wave simulations. Additionally, the long-wave theory is strictly valid for ≡ / ≪ 1 but we solved the evolution equation in a finite domain with a small but finite value of . We have not solved the full axisymmetric Stokes problem here, which could be used to validate the long-wave model.
Our model predicts that viscoplastic rheology can significantly alter the evolution of a layer coating a cylindrical tube. When the layer is thin, key aspects of the dynamics and the final shape of the layer can be understood by studying marginally-yielded static solutions. When the layer has finite thickness, the critical thickness required to form a liquid plug can depend strongly on the capillary Bingham number. These results have implications for modelling real-world problems where the coating liquid has a yield stress, such as models of airway closure, particularly in the context of diseases which alter mucus rheology.
A.2. Energy in the long-wave system In the rest of Appendix A, we use subscripts to denote derivatives. The non-dimensionalised expression for energy, , in the long-wave system is given in (2.21). Differentiating (2.21) with respect toˆ , then inserting the evolution equation (2.16), integrating by parts, and using the boundary conditions (2.18), gives whereˆ is defined in (2.14). Expandingˆ using the definition in (2.16) gives where the functions 1 ( , Ψ) and 2 ( , Ψ) are defined in (2.17). From (A 7), ˆ 0 if 1 ( , Ψ) 0 and 2 ( , Ψ) 0 for all ( , Ψ) ∈ D ≡ {( , Ψ) : 0 Ψ 1}. To prove that this is the case, first note that 1 and 2 are both monotonically decreasing in , for ( , Ψ) ∈ D. This can be seen from noting that The functions 1 (Ψ) and 2 (Ψ) are non-negative for all 0 Ψ 1 (figure 9a) so the derivatives (A 8) are both non-positive for all ( , Ψ) ∈ D. This implies that if 1 and 2 are non-negative on the boundary Ψ = of D, then they are non-negative everywhere in D. Setting Ψ = in (2.17), we find the functions The functions 1 ( , ) and 2 ( , ) defined in (A 10). All four functions are non-negative, which is used to prove that ˆ 0.
A.3. Energy in the thin-film system We can also show directly from the thin-film equations that energy is non-increasing. The interfacial energy in the thin-film limit is (2.31), which when differentiated with respect to gives where we have used integration by parts, the boundary conditions (2.28) and the evolution equation (2.27). Noting that 0 , (A 11) immediately implies 0.

Appendix B. Derivation of the long-wave evolution equation
Starting from the governing equations and boundary conditions in the Stokes system (2.1)-(2.8), we derive the long-wave evolution equation (2.13)-(2.18). First, we rewrite (2.1)-(2.8) in terms of the non-dimensionalised and scaled variables (2.11) and (2.12). The Stokes equations (2.2) become The non-zero components of the shear-rate (2.1) become and the second invariants of stress and shear-rate are The constitutive relation (2.3) becomes whereˆ ≡ / . From (2.4), the wall boundary conditions arē and from (2.5)-(2.7), the free-surface boundary conditions are The symmetry boundary conditions (2.8) become A description of the flow where the fluid is yielded is first derived, then regions of unyielded fluid can subsequently be identified. Until otherwise stated, we assume >ˆ . Separate expansions are made for the relevant variables in the shear-dominated region, 1, and in the plug-like region, < . Quantities in the shear-dominated region are denoted by the superscript (·) , and quantities in the plug-like region by (·) . In the shear-dominated region, let = 0 + 1 + . . . ,¯ =¯ 0 + ¯ 1 + . . . , = 0 + 1 + . . . , The leading-order second invariant of stress is then 0 = | 0 |. After truncating at leading order, the horizontal momentum equation (B 1c) becomes 0 = − ¯ ¯ + 1 0 in 1, (B 10) and the vertical momentum equation (B 1b) implies¯ =¯ (¯ ,¯ ) in 1. The sheardominated solution is only valid where | 0 | >ˆ , so we identify that | 0 | = | 0 | =ˆ at = . After integrating (B 10) in , we enforce | 0 | =ˆ at = to get The shear component of the constitutive relation (B 4), at leading-order, gives In both (B 11) and the second equality of (B 12), we have used the fact that the shear stress, 0 , must have the same sign as the pressure gradient. Combining (B 11) and (B 12) gives Integrating (B 13), and enforcing the no-slip condition at = 1, finally gives in 1. In the plug-like region, < , we make an expansion of the same form as (B 9), = 0 + 1 + . . . ,¯ =¯ 0 + ¯ 1 + . . . , = 0 + 1 + . . . , but we assume the leading-order axial velocity is independent of , so 0 = 0 (¯ ,¯ ). This means that = ( ), and so = ( ). Thus, the shear component of the constitutive relation (B 4) implies with (B 3b) giving the leading-order second invariant, After truncating, the horizontal momentum equation (B 1c) becomes The expression (B 21) is valid up to the point at whichˆ = | 0 |, which must coincide with = . Hence, using (B 20), we arrive at the definition, Note that from (B 22), always holds. The axial velocity in the shear-dominated region (B 14) is matched to 0 by equating them at = , which gives Then (B 14) and (B 23) provide the complete expression for axial velocity, 0 , across the whole layer, 1, in regions where the fluid is yielded. Using the same argument as Balmforth & Craster (1999), we identify that if (¯ ,¯ ) 1 then the shear-dominated region does not exist, and the boundary conditions (B 5) imply the plug-like region must be stationary, 0 = 1 = 0, so the fluid is unyielded. If we replace with Ψ(¯ ,¯ ) ≡ min(1, ) in (B 14) and (B 23), then 0 = 0 for < Ψ and 0 = 0 for Ψ 1 hold both where the fluid is yielded and where it is unyielded. The leading-order axial flux is then which when evaluated gives the expression in (2.16). Finally, the kinematic boundary condition (B 6a) and mass conservation (B 1a) are combined to give the evolution equation (2.16). When we present the evolution equation and relevant definitions in (2.13)-(2.18), we write them in terms of the unscaled variables ,ˆ ,ˆ ,ˆ , , i.e. we view the system in a frame unscaled by the small aspect ratio . However, the limit in which the theory is formally valid remains ≪ 1.
the vector = ( 1 , 2 ) is the solution to L = 0 with boundary conditions B = 0. Figure  10(a) shows the computed solution. Note that the amplitude of 1 is free, so to solve we choose an arbitrary value by setting 1 ( ) = 1. We now note the linear ODE problem defined by (C 3) has an associated adjoint problem. The adjoint operator, L † , and boundary operator, B † , are defined via the relation † , L = L † † , , where the inner product is defined as , = ∫ 0 d for vectors , . This The adjoint solution, † = ( † 1 , † 2 ) , satisfies L † † = 0 with boundary conditions B † † = 0. The function † 2 is constant; its value is arbitrary but sets the amplitude of † 1 . Figure 10(a) shows the computed solution for † 1 where we set † 2 = 1. To find A ( ), first note that the evolution equation (2.27) implies H 1 = ( 4 ). We introduce a slow timescale = 3 and let A = A ( ). An evolution equation for A ( ) is found by solving the ( 2 ) problem. At ( 2 ), we get which is the evolution equation for A ( ). Notice that (C 9) is independent of the amplitude of † 1 , so the value of † 2 is truly arbitrary. After computing and † , the coefficients in (C 9) are found using numerical integration. Equation (C 9) has two fixed points at A ≈ ±2.20. Figure 10(b) is a solution of (C 9) with initial conditions A (0) = −2.19, showing that the solution evolves away from the negative fixed point towards the positive one, suggesting the former is unstable and the latter stable. Fixed points in A correspond to static solutions for . Since 1 ( ) > 0 (figure 10a), the negative fixed point must correspond to the lower-branch static solution in figure 4(a), since then H 1 < 0 so max < max 0 ( ; * ). Similarly, the positive fixed point must correspond to the upper-branch solution since H 1 > 0. This confirms that, at least near = * , the lower-branch solutions are unstable and the upper-branch solutions are stable. The Region I problem is equation (D 4) subject to boundary conditions (D 7a) and (D 15). When solving the Region I problem, we define a small constantˆ and solve (D 4) in the domain 0 − −ˆ . We enforce (D 15) by settingĥ 0 ( − −ˆ ) = 8/3ˆ 3/2 and ℎ ′ 0 ( − −ˆ ) = √ 6ˆ 1/2 . We chooseˆ sufficiently small that the solution in the rest of the domain is insensitive to its exact value.
When solving the Region II problem to findh 0 , we define a large constant ∞ and solve (D 5a) in the finite domain − ∞ < < ∞ . Informed by (D 9a) and (D 14), we enforce the boundary conditions,h 0, We choose ∞ large enough that the solution is insensitive to its exact value away from the boundaries = ± ∞ . Figure 11(a) shows the composite solution, which we call ( ; ), for = 10 −3 . To compute this solution, we usedˆ = 10 −4 and ∞ = 600. The value 1 ≈ −0.105 computed with this solution completes the Region III solution (D 13). The constant 1 depends on , which here is = √ 2 . For clarity, the solutions displayed in figure 11 are truncated shortly after the points where they overlap. This also means that the parts of the solutions displayed are away from the boundaries so entirely independent of the values ofˆ and ∞ chosen. Figure 11(b) shows the capillary pressure, ≡ − − , , of the composite solution.

Appendix E. Static solutions of the long-wave evolution equation
Following our approach in §3.3, we look for marginally-yielded static solutions, = 0 ( ; , ), of the long-wave equations (2.13)-(2.18). The static shapes 0 ( ; , ) are solutions to = 1 where is defined in (2.13). As in the thin-film analysis, we assume that the pressure is monotonic, soˆ < 0. The ODE = 1 can be rearranged to givê (1 − 2 0 ) = −2 2 , (E 1) whereˆ is defined in (2.14) and 2 =ˆ . We solve (E 1) subject to boundary conditions 0, (0; , ) = 0, ( ; , ) = 0, and the volume conservation condition, The problem is solved using a boundary value problem solver in M . We define the layer thickness, 0 ( ; , ) ≡ (1 − 0 )/ , to aid discussion and comparison with the thin-film static solutions computed in §3.3. Figure 12 shows solutions for = 0.1, 0.12, 0.14. Figure 12(a) shows that, like in the thinfilm case, we find an upper and a lower branch of solutions for each , and a bifurcation point = * such that no solutions exist for > * . Note that the location of the bifurcation now depends on . The boundary value problem solver is generally able to compute the whole lower branch and most of the upper branch of solutions, except for very small . The upper-branch solutions are very singular for small , with an increasingly large jump in 0, around the minimum in 0 , which makes computing them difficult. Figure 12(b) shows plots of 1 − 0 (0; , ). In §3.4, we show that the same quantity from the thin-film static solutions has particular significance in determining the outcome of an evolving thin layer. Figure 8 shows that it has much less physical significance in the long-wave problem when 0.12. We argue that this is because the evolving layer generally does not select a static shape which is a solution to (E 1). Instead, either a plug forms, or the layer rigidifies near = 0 early in the evolution which leads to a different static two-collar solution being selected.