Surfactant amplifies yield-stress effects in the capillary instability of a film coating a tube

To assess how the presence of surfactant in lung airways alters the flow of mucus that leads to plug formation and airway closure, we investigate the effect of insoluble surfactant on the instability of a viscoplastic liquid coating the interior of a cylindrical tube. Evolution equations for the layer thickness using thin-film and long-wave approximations are derived that incorporate yield-stress effects and capillary and Marangoni forces. Using numerical simulations and asymptotic analysis of the thin-film system, we quantify how the presence of surfactant slows growth of the Rayleigh-Plateau instability, increases the size of initial perturbation required to trigger instability and decreases the final peak height of the layer. When the surfactant strength is large, the thin-film dynamics coincide with the dynamics of a surfactant-free layer but with time slowed by a factor of four and the capillary Bingham number, a parameter proportional to the yield stress, exactly doubled. By solving the long-wave equations numerically, we quantify how increasing surfactant strength can increase the critical layer thickness for plug formation to occur and delay plugging. The previously established effect of the yield stress in suppressing plug formation [Shemilt et al., J. Fluid Mech., 2022, vol. 944, A22] is shown to be amplified by introducing surfactant. We discuss the implications of these results for understanding the impact of surfactant deficiency and increased mucus yield stress in obstructive lung diseases.


Introduction
Pulmonary surfactant plays a crucial role in the healthy function of the lungs.By lowering the surface tension at the interface between air and the liquid that lines lung airways, surfactant helps to prevent unwanted collapse of small airways and alveoli during breathing (Milad & Morissette 2021).Surfactant deficiency is likely to contribute to the increased prevalence of airway obstructions in diseases such as asthma (Hohlfield 2002) and cystic fibrosis (Tiddens et al. 2010).Mucus, the main component of the airway surface liquid, is a complex fluid exhibiting properties such as viscoelasticity, shear-thinning and a yield stress (Hill et al. 2022).In various obstructive diseases, and particularly in cystic fibrosis, the mucus typically has altered rheology, including a significantly increased yield stress compared with mucus in healthy lungs (Patarin et al. 2020).This provides motivation for this study into the effect of insoluble surfactant on the instability of a viscoplastic liquid coating the interior of a cylindrical tube, which is a simple model for the flow of mucus that can lead to airway closure.Additional applications can be found in related interfacial flows of viscoplastic fluids from engineering and industry, where surfactants may also be present (Mitsoulis 2007;Glasser et al. 2019;Ahmadikhamsi et al. 2020).
The instability of a liquid film coating the interior of a cylindrical tube has been widely studied in the case where the liquid is Newtonian.When the volume of fluid in the layer is small, annular liquid collars will form on the tube wall, and when the volume is large enough, liquid plugs form in the tube (Everett & Haynes 1972).Hammond (1983) derived an evolution equation for the motion of a thin layer, and presented numerical and late-time asymptotic solutions showing annular collars forming and fluid slowly draining out of thin regions between them.This thin-film theory was then extended to describe the motion of thick films by Gauglitz & Radke (1988), who deduced numerically that an approximate minimum thickness of 12 % of the tube radius is required for plug formation to occur.Otis et al. (1990) were the first to extend the theory to include the effect of insoluble surfactant at the air-liquid interface.They found that growth of the instability and plug formation is delayed by the surfactant.Moreover, if the surfactant is strong, then Marangoni forces effectively immobilise the interface, increasing the time scale for the evolution of the layer by a factor of four compared with when there is no surfactant.This factor of four decrease in the growth rate had been previously identified by Carroll & Lucassen (1974) in the related instability of a thin liquid layer coating the exterior of a cylindrical filament.Results from the thick-film model of Ogrosky (2021) suggest that, whilst this factor is very close to four for layers with thicknesses close to the critical value for plug formation, it may be increased for thicker layers.Halpern & Grotberg (1993) investigated the effect of surfactant on the evolution of a layer coating the interior of an elastic tube, and also found slowing of the dynamics and a delay to plug formation, except when the tube stiffness was very low, in which case the impact of including surfactant was minimal.They argued that this slowing implies an increase in the critical thickness required for plug formation, since simulations were run to a fixed finite time.This observation has relevance to mucus plug formation in airways, which typically form within the time scale of a single breath cycle.Experimental results have confirmed the decreased growth rates and increased times for plug formation to occur due to surfactant (Cassidy et al. 1999).Computational fluid dynamics (CFD) simulations have also been conducted which, unlike quasi-one-dimensional models, can capture the post-coalescence phase of plug formation as well as the pre-coalescence phase (Romanò, Muradoglu & Grotberg 2022).It was found, again, that introducing surfactant delayed plug formation, and also that it decreased the stress on the tube wall during plugging.The shear stress exerted on the tube wall is a physiologically important variable in airway closure models as it has been shown that epithelial cell damage can be caused by sufficiently large stresses being exerted on the airway wall (Huh et al. 2007).
There has been some attention on the effect of non-Newtonian liquid rheology on this flow in the case where there is no surfactant present.Halpern, Fujioka & Grotberg (2010) studied the effect of viscoelasticity, showing that the time for a plug to form can be decreased by increasing the Weissenberg number if the layer is sufficiently thick.Romanò et al. (2021) also investigated a viscoelastic version of the flow, using CFD, and found that elastic effects can induce a significant peak in the shear stress on the tube wall in the post-coalescence phase of plug formation.Erken et al. (2023) took a similar approach but with an elastoviscoplastic model for the liquid layer, showing that elastic effects can also impact how the fluid yields, particularly around and immediately after plug formation.Shemilt et al. (2022) derived reduced-order models to study the effect of viscoplastic rheology on the dynamics of thin films and of thick films in the lead-up to plug formation.It was found that increasing the capillary Bingham number, B (a parameter proportional to the liquid yield stress), can suppress instability by rigidification of the layer, reduce deformation when there is instability and significantly increase the critical layer thickness required for plug formation to occur.A model describing the evolution of thicker layers was developed using a long-wave approximation, and a simplified evolution equation was deduced using thin-film theory.The viscoplastic long-wave theory used is closely related to the planar thin-film theory exposed by Balmforth & Craster (1999).The structure of the flow is qualitatively the same in long-wave and thin-film theories: where the shear stress exceeds the yield stress, the fluid is fully yielded and the flow is shear-dominated, but where the shear stress is less than the yield stress, the flow is plug-like with no shear flow at leading order.However, due to changes in the surrounding flow, some regions with plug-like flow must still deform.In these regions, the yield stress is exceeded by an asymptotically small amount (Balmforth & Craster 1999) and the regions are referred to as 'pseudo-plugs' (Walton & Bittleston 1991).Whilst the long-wave theory is similar to thin-film theory, it differs in a few crucial ways: the layer is not assumed to be thin relative to the tube radius, additional terms are retained in the evolution equations that more accurately capture the curvature of the geometry, and the exact expression for the curvature of the interface is used instead of the linearised version used in the thin-film theory.Thus, long-wave theory provides a composite approximation to the dynamics of a thick film, which is accurate where the layer is thin and also describes well the regions of the flow that are approximately capillary static.Equivalent, or very similar, approximations have been used previously to study this flow when the liquid layer is viscoelastic (Halpern et al. 2010) or Newtonian (Gauglitz & Radke 1988;Johnson et al. 1991;Camassa, Ogrosky & Olander 2014, 2017), including when insoluble surfactant is present (Otis et al. 1990(Otis et al. , 1993;;Ogrosky 2021).
Viscoplastic thin-film theory has been used to describe various other interfacial flows where surface tension plays a key role.Balmforth, Ghadge & Myers (2007) studied the surface-tension-driven fingering instability in flow down an inclined plane, showing that the yield stress has a stabilising effect.Jalaal & Balmforth (2016) developed a model using thin-film theory for a propagating bubble through a tube filled with viscoplastic fluid, which they validated against CFD simulations.Increasing the Bingham number was found to increase the thickness of the film between the bubble and the tube wall.Thin-film theory compared well with simulations for low Bingham numbers but less well when the film thickness increased.Thin-film models for the axisymmetric spreading of droplets (Jalaal, Stoeber & Balmforth 2021) and the spreading of long extruded filaments (van der Kolk, Tieman & Jalaal 2023) have been developed and used to predict the distances reached by spreading fronts, which decrease as the capillary Bingham number is increased.Whilst these studies all addressed the effect of viscoplastic rheology on capillary phenomena, surfactant effects were not incorporated in any of the models.Craster & Matar (2000) used thin-film theory to model surfactant-driven flow in viscoplastic films, demonstrating that after Marangoni forces cause a spreading front to develop, the yield stress can rigidify the layer before it returns to a uniform height profile.This study did not, however, include the effects of capillary forces.To the authors' knowledge, viscoplastic thin-film theory has not previously been used to study any flows where both capillary and Marangoni forces are present.
In this study, we develop a model for the evolution of a liquid film coating the interior of a tube, where the flow is driven by surface tension but is also influenced by Marangoni forces.Our aim is to quantify the effects of surfactant on the capillary instability and on the previously established effects of the liquid's yield stress (Shemilt et al. 2022).We will derive evolution equations for the layer height and surfactant concentration using long-wave theory, and subsequently deduce a simpler version of these equations that is valid in the thin-film limit.To focus attention on the interaction between Marangoni, capillary and yield-stress effects, other phenomena relevant to airway modelling are not included in the model.The Bingham constitutive model is used, which does not include rheological properties such as shear-thinning, viscoelasticity or thixotropy, but allows us to focus attention on yield-stress effects.Surface tension is assumed to vary linearly with surfactant concentration.For pulmonary surfactants, this relation is nonlinear (Schürch, Bachofen & Possmayer 2001), but the linear model captures key Marangoni effects while being amenable to detailed analysis.This choice is in keeping with previous models (Halpern & Grotberg 1993;Ogrosky 2021).Moreover, gravity is assumed negligible, the tube wall is assumed rigid and the air in the centre of the tube is assumed passive.Numerical solutions of both the long-wave and thin-film equations will be used to elucidate the new features of the dynamic evolution that arise due to the presence of surfactant, and also to systematically explore parameter space.We will exploit the relative simplicity of the thin-film equations to study the behaviour of the layer in a late-time limit and in the limit of large Marangoni number (when surfactant is strong).By computing and analysing solutions of the long-wave equations, we will quantify how surfactant alters the dynamics leading to plug formation and the critical thickness required for plugging to occur.In particular, we will show how surfactant acts synergistically with the yield stress to stabilise the liquid layer.
The paper will proceed as follows.In § 2, we derive evolution equations for the layer thickness and surfactant concentration, with the long-wave equations given in § 2.3 and the thin-film equations in § 2.4.Solution methods will then be briefly discussed in § 2.5.Results from the thin-film system will be presented in § 3.An example numerical simulation will be discussed in § 3.1, late-time asymptotic analysis of the thin-film equations will be presented in § 3.2, and the effect of varying the Marangoni number on the stability and evolution of the layer will be systematically addressed in § 3.3.Results for thick films from the long-wave theory will be given in § 4.An example numerical simulation will be examined in § 4.1, a discussion of the behaviour when the Marangoni number is large will be given in § 4.2 and the effect of surfactant on the critical thickness required for plug formation will be explored in § 4.3.Finally, in § 5, there will be a discussion of the significance of the results, particularly in relation to modelling airway closure in the lungs.

Governing equations and boundary conditions
We consider a rigid circular cylindrical tube lined on its interior by a layer of viscoplastic liquid, with a gas in the centre of the tube and insoluble surfactant present at the gas-liquid interface (see figure 1).We assume the system is axisymmetric so it can be described using cylindrical coordinates (r * , z * ).The tube has radius a and the interface is located at r * = R * (z * , t * ).The liquid layer has velocity u * (r * , z * , t * ) = u * r + w * ẑ and pressure p * (r * , z * , t * ) relative to the gas pressure, which is assumed to be spatially uniform.
The liquid is incompressible and we assume inertia is negligible, so conservation of mass and momentum imply where τ * is the deviatoric stress tensor.The boundary conditions at the interface are the kinematic condition, and the stress condition, where σ * is the surface tension, κ * is the curvature of the interface, n is the unit normal to the interface directed away from the liquid layer and ∇ * s is the surface gradient operator.At the tube wall, there is no slip and no penetration, (2.4) The liquid is assumed to be a Bingham fluid with constitutive relation, where η is a viscosity, τ y is the yield stress, γ * = ∇ * u * + ∇ * u * T is the shear-rate tensor, and τ * and γ * are the second invariants of τ * and γ * , respectively.The second invariant of a tensor T is defined as T = 1 2 T ij T ij .We take the surface tension of the interface to be a linear function of the surfactant concentration, where σ 0 and Γ 0 represent the constant values of the surface tension and surfactant concentration in an unperturbed state, and K > 0 is a constant.We assume that the surfactant is insoluble so its motion is governed by the transport equation (see, e.g.Stone 1990) where u * s is the velocity along the interface.In (2.7), we have neglected diffusion of surfactant to focus on the limit of advection-dominated transport.In lung airways, although there is significant variation in measurements of the properties of surfactant and mucus, we typically expect surfactant diffusivity to be small and that advection will be the dominant transport mechanism (Craster & Matar 2000;Lai et al. 2009;Chen et al. 2019).At the lateral boundaries of the domain, we impose symmetry boundary conditions, which enforce that the first derivative of the layer height is zero and that there is no flux of fluid or surfactant across the side boundaries.

Surfactant amplifies yield-stress effects
The constitutive relation (2.5) is is a capillary Bingham (or plastocapillarity) number (Jalaal et al. 2021;van der Kolk et al. 2023).The equation of state (2.6) becomes where is a Marangoni number.The transport equation (2.7) can be expressed in conservative form (Halpern & Frenkel 2003) as where w s is the axial component of the surface velocity.We now examine simplified versions of (2.10)-(2.21),first in a long-wave limit, appropriate for thicker films, and then in a more restrictive thin-film limit.
2.3.Long-wave theory We consider the system (2.10)-(2.21) in a long-wave limit by defining a characteristic axial length scale, L = a/δ, where δ 1.We define stretched variables, where the choice of scalings for z and L arise from the long-wave approximation, the scaling for τ is then chosen so that the radial gradient of shear stress balances the axial pressure gradient in (2.12), the scaling for w is chosen so that the shear stress balances with the corresponding term in the strain-rate tensor in (2.17), u is scaled such that the mass conservation equation (2.10) balances and, finally, the scaling for t allows all terms in the kinematic boundary condition (2.14) to balance.(The scalings (2.22) correct those printed in (2.11) of Shemilt et al. (2022) but the equations derived there are still correct and consistent with what we derive here.)We then truncate the governing equations (2.10)-(2.21) at leading order in δ, with the exception of (2.16) where we retain the exact curvature, which is a commonly used device to improve accuracy in near-static regions of the flow (Gauglitz & Radke 1988;Halpern & Grotberg 1992, 1993;Halpern et al. 2010;Ogrosky 2021;Shemilt et al. 2022).The leading-order governing equations and boundary conditions are given in Appendix A, where we also detail the derivation of the long-wave evolution equations, which are presented in the remainder of this section.As has been done previously in similar problems (Camassa et al. 2012;Camassa & Ogrosky 2015;Ogrosky 2021;Shemilt et al. 2022), we will present the long-wave equations here in terms of the unscaled variables defined in (2.9) instead of the scaled variables (2.22), but the equations still represent the leading-order theory in the limit δ 1.To understand the structure of the evolution equations, it is important to first understand the structure of the flow and the ways in which the layer can yield.Where the magnitude of the shear stress is larger than the yield stress, |τ rz | > B for some R(z, t) < r < 1, the fluid is fully yielded and the flow is shear-dominated.Where |τ rz | ≤ B for some R(z, t) < r < 1, the flow is plug-like and the leading-order axial velocity is independent of r, i.e. w = w p (z, t).If w p is also independent of z in a plug-like region, then it is a rigid plug.If not, then that plug-like region is deforming axially so it must be yielded.In those yielded plug-like regions, the normal stresses become leading-order in δ and combine with the shear stress such that the yield condition is just met; such a region is referred to as a 'pseudo-plug' (Walton & Bittleston 1991).This general structure is the same as in viscoplastic thin-film flows, as detailed by (Balmforth & Craster 1999).For the remainder of this section, subscripts will be used to denote derivatives.
Where the fully yielded regions, pseudo-plugs and rigid plugs occur in the liquid film depends on the competition between surface and bulk forcing, via M Γ z and p z , and their sizes relative to B. In Appendix A, we show that, in the long-wave limit, the capillary pressure is given by (2.23) where the curvature of the interface, κ, is defined in (2.16), and that the shear stress is given by (2.24) By noting how τ rz depends on r in (2.24), it can be deduced that there are five qualitatively different types of yielding that can occur in the layer, corresponding to the five possible combinations of fully yielded regions, rigid plugs and pseudo-plugs that can exist in R ≤ r ≤ 1.We define two internal surfaces, r = Ψ − (z, t) and r = Ψ + (z, t), which separate fully yielded and plug-like regions.The fully yielded regions (where |τ rz | > B) occupy R ≤ r ≤ Ψ − and Ψ + ≤ r ≤ 1, and the plug-like region (where Note that each of these three regions may not always exist, in which case the width of the region will be zero.The five possible types of yielding are as follows and are illustrated in figure 2(a).
(I) Internal pseudo-plug (R < Ψ − < Ψ + < 1).Here, fully yielded regions adjacent to the wall and the interface are separated by a pseudo-plug.(II) Surface pseudo-plug (Ψ − = R < Ψ + < 1).Here, the pseudo-plug extends to the interface and the only fully yielded region is adjacent to the wall.(III) Fully yielded (Ψ − = Ψ + = R or 1 = Ψ − = Ψ + ).In this case, there is no plug-like region and the whole layer is fully yielded.(IV) Near-wall plug (R < Ψ − < 1 = Ψ + ).Here, there is a rigid plug adjacent to the wall and the only fully yielded region is adjacent to the interface.(V) Fully rigid (Ψ − = R and Ψ + = 1).Neither yielded region exists and the whole layer is a rigid plug.
Hewitt & Balmforth (2012) identified these same five yielding types in the flow of a thin film between two moving solid surfaces.In that problem, there are simple criteria to determine which yielding type occurs based on the size of the shear stress at the two solid boundaries.Here, the additional complexity of the long-wave theory compared with the thin-film theory means there are not such simple criteria.Instead, we derive the following The fives types of yielding that can occur in the layer.Plug-like regions are shown in grey and fully yielded regions in white.Typical axial velocity profiles are also sketched.Below are maps of parameter space showing where these yielding types occur in (b) the long-wave system and (c) the thin-film system.When plotting the map in panel (b), we treat R as a fixed parameter to focus on variation with p z and M Γ z .In panel (b), as in § 2.3, the parameter map is plotted in terms of the unscaled variables (2.9).In panel (c), the map is plotted in terms of the scaled thin-film variables introduced in (2.37).Along the dashed line in panel (c), the surface velocity is exactly zero, ws = 0.
expressions for Ψ ± , from which the type of yielding can be deduced.We find where, if p z / = 0, we have Given values of R, B, p z and M Γ z , the type of yielding can then be deduced from (2.26).A more intuitive illustration of when the yielding types I-V occur is given by figure 2(b).It shows how the yielding type depends on the size of the capillary stress relative to the yield stress, via p z / B, and on the size of the Marangoni stress relative to the yield stress, via M Γ z / B. We plot figure 2(b) assuming R is fixed, even though, in general, it will vary with z.We now briefly survey ( p z / B, M Γ z / B)-space.Starting in the upper left region of figure 2(b), where capillary and Marangoni stresses are both large but with opposite signs, there is yielding of type I. Suppose M Γ z / B is then decreased, so that we cross the line M Γ z / B = 1.Then the shear stress at the interface has dropped below B so there can no longer be a fully yielded region at the interface, and the layer then exhibits yielding of type II.Decreasing M Γ z / B further so that M Γ z / B < −1, the Marangoni stress at the interface again exceeds the yield stress so there must be yielding at the interface, but it now acts in the same direction as the capillary stress.The layer then exhibits yielding of type III.Now increasing p z / B, we remain in yielding type III until we reach the line After crossing this line, the capillary stress is no longer strong enough to yield the fluid adjacent to the wall, so a rigid plug develops there and we have yielding of type IV.Increasing p z / B yet further, so that we cross the line is large enough that the lower fully yielded region appears again so we have returned to yielding of type I, but with the flow in the opposite direction to when we started.Symmetry means that if we proceed in the same fashion around the other half of the plane, the yielding transitions will be the same as just described.There are two regions in figure 2 that we have not yet discussed.When p z / B and M Γ z / B are both small, there is no yielding, so there is a region with yielding of type V around the origin.Finally, yielding of type I can also be observed in two small regions near ).The existence of these relies on the shear stress in the long-wave theory being nonlinear, so they do not exist in the simpler thin-film limit (figure 2c).
Given the expressions for Ψ ± in (2.25) and (2.26), we can now present the long-wave evolution equations, which are derived in Appendix A. The evolution equation for the interface position, R, is where the axial volume flux is given by (2.28) with ) If M = 0, or equivalently if there is no surfactant, Γ = 0, then Ψ − = R and (2.27)-(2.29)reduce to the evolution equation for the surfactant-free problem (correcting a typographical sign error in the expression for Q given in (2.16) of Shemilt et al. 2022).
The surfactant transport equation is where the surface velocity is with Equations (2.27)-(2.32)are solved subject to the boundary conditions which completes the long-wave system of equations.Finally, by evaluating the shear stress (A7) at r = 1, we can deduce an expression for the stress exerted on the tube wall, (2.34) 2.4.Thin-film theory We now consider the system in a thin-film limit to derive simpler evolution equations that are more amenable to detailed analysis.In the thin-film theory, we assume that |1 − R| 1, but no longer require that δ 1.Rather than presenting a derivation of the thin-film equations from (2.10)-(2.21),here we will show how they can be deduced directly from the long-wave equations (2.27)-(2.33).The flow structure is qualitatively the same and the same five possible yielding types I-V also occur.
We assume that where 1 and H is the scaled layer thickness.Similarly, we define the boundaries between fully yielded and plug-like regions, (2.36) Since Ψ + ≥ Ψ − , the definitions (2.36) mean Y + ≥ Y − .Other relevant variables are rescaled by defining where assuming for now that p z / = 0.As in § 2.3, we can present criteria for each of the five yielding types to occur based on the values of Y ± (as in figure 2a): I, internal pseudo-plug (0 The evolution equation (2.27), to leading order in , becomes where the scaled axial volume flux is As can be seen in figure 2(c), when pz = 0, the only type of yielding possible is type III, which occurs if The surfactant transport equation (2.30), at leading order, becomes where the scaled surface velocity is

Surfactant amplifies yield-stress effects
This line is included in figure 2(c) and will form part of our later discussion.The lateral boundary conditions (2.33) reduce in the thin-film limit to (2.45) Equations (2.40)-(2.44)with boundary conditions (2.45) comprise the thin-film system.The wall shear stress (2.34), in the thin-film limit, becomes (2.46) 2.5.Solution methods When solving the thin-film or long-wave equations numerically, we use initial conditions with a uniform surfactant concentration across the layer and a perturbation to the interface height with wavelength 2L.For the thin-film equations, these initial conditions are and for the long-wave equations, the equivalent conditions are where A is the amplitude of the initial perturbation and is the ratio of the average layer thickness to the tube radius when A = 0.The constant term in (2.48) ensures that the volume of fluid is independent of A. As noted for the surfactant-free problem (Shemilt et al. 2022), when the fluid is viscoplastic, there is no linear instability since a finite-amplitude initial perturbation is required to yield.
For ease of comparison with the surfactant-free problem (Shemilt et al. 2022) and other related studies (Gauglitz & Radke 1988;Halpern et al. 2010), we solve the equations in a domain of length L = √ 2π.This is the domain length that can accommodate the most unstable mode in the Newtonian linear stability analysis (Hammond 1983).The initial perturbation applied to the layer then corresponds to the one unstable Fourier mode that exists in the domain.Although this value of L is not necessarily of unique importance in the viscoplastic problem, we do not observe any qualitative changes in the results when L is changed by relatively small amounts.As in the surfactant-free problem, since there is nothing in the model equations that could break symmetry around the lateral boundaries of the domain, the boundary conditions (2.33) and (2.45) yield the same solutions as would be found using periodic boundary conditions.Symmetry rather than periodic boundary conditions allows a finer grid to be used for the spatial finite differencing at the same computational expense, since the computational domain is shorter.
To solve both the long-wave equations and the thin-film equations numerically, we use a regularisation introduced by Jalaal (2016), which was used for the surfactant-free problem (Shemilt et al. 2022) and has also been used for several other viscoplastic thin-film problems (Balmforth et al. 2007;Jalaal et al. 2021).We define Ŷ± = max(Y min , Y ± ) and Ψ± = min(1 − Y min , Ψ ± ) where Y min is a small parameter, and replace Y ± and Ψ ± with Ŷ± and Ψ± , respectively, in the equations.We choose Y min = 10 −8 for all simulations, which is small enough that the exact value does not affect the results.To solve the resulting equations, the domain 0 ≤ z ≤ L is discretised into a grid of N evenly spaced points (we typically use N = 200), the spatial derivatives are approximated by second-order central finite differences and then the equations are time-stepped using an ordinary differential equation (ODE) solver in MATLAB.

Time evolution of a surfactant-laden thin film
Figure 3(a) shows snapshots from a sample numerical solution of the thin-film system.At t = 0, the initial perturbation amplitude, A, is large enough that the layer is yielded in a region in the centre of the domain with a fully yielded region at the base of the layer.At very early times, this yielded region spreads to cover the whole domain by t = 5.In this very-early-time period, Γ z remains small, so the dynamics are essentially uninfluenced by the presence of surfactant.The behaviour is qualitatively the same as in the early-time yielding period described previously for the surfactant-free problem (Shemilt et al. 2022).Figure 3(b) shows that there is a delay in the initial growth of the instability compared with a Newtonian (B = 0) simulation, and that the delay is approximately equal in the surfactant-free viscoplastic simulation.After this early-time period, however, figure 3(b) shows that the growth rate of the instability is reduced by the presence of surfactant.
Between t = 5 and t = 80 (figure 3a), significant gradients in surfactant concentration develop, and at t ≈ 80, a fully yielded region at the interface, where MΓ z exceeds B, appears near the centre of the domain.This yielded region grows and has extended across the whole interface by t = 500.During the intermediate period, 80 < t < 500, the two lateral edges of the upper yielded region propagate towards the side boundaries; these coincide with the location of two gradually steepening travelling wave fronts in Γ z .At t = 130, the right-travelling wave in Γ z near z = 4 resembles a discontinuous shock wave.The sudden decrease in Γ z across the shock results in a rise in Y − just ahead of the shock where the Marangoni force is weaker.Similarly, at t = 250, a shock-like discontinuity has developed in Γ z in the left-travelling wave as it approaches z = 0, and a small rise in Y − can be observed ahead of the wave.
Once we observe that a shock has developed in Γ z , we can use (2.43) to derive a Rankine-Hugoniot condition (see Appendix B), which provides the relation for the shock propagation speed, u s , in terms of the sizes of jumps in Γ z and ( ws Γ ) z across the discontinuity.We have already noted how the jump in Γ z across the shock affects the yielding behaviour on either side, but (3.1) also shows that the yielding behaviour, via ws , influences the shock propagation, highlighting the coupling between rheology and surfactant transport.Although our numerical method does not actively track the shock location, the speed of shock propagation observed in the simulations agrees well with u s calculated via (3.1) (data not shown here), providing evidence that the numerics accurately capture the behaviour around the shock.Throughout the intermediate-time period described above, the regions ahead of the shock waves exhibit yield type II (surface pseudo-plug) while the central region behind the shocks exhibits yield type I (internal pseudo-plug).At t = 80, t = 100 and t = 130, at both side boundaries, the solution in (H pz /B,MΓ z /B)-space approaches (H pz /B, MΓ z /B) = (0, −1).This indicates that near the side boundaries, capillary forces are small and Marangoni forces dominate.Unlike in the clean problem, where Y − → 0 as z → 0 or z → L, here, due to the presence of Marangoni forces, Y − is not necessarily zero at the side boundaries.The value of Y − in the limits z → 0 or z → L during this period is given by lim z→0,L (H − MΓ zz /p zz ), which can be seen to take a positive value in figure 3 2000), the whole solution is close to the line MΓ z + 1 2 H pz = 0.This indicates that H − Y + ≈ Y − , and so ws ≈ 0 from (2.44), meaning the interface is approximately immobilised at late times.As the layer approaches its final static shape, all points in the domain converge towards (H pz /B, MΓ z /B) = (−2, 1) (figure 3a, t = 2000).Hence, when the layer reaches its final static shape, the capillary stress is twice as strong as the Marangoni stress and they act in opposite directions, with the resultant magnitude of stress exactly equal to B. From this, we can deduce that at late times, the layer approaches a marginally yielded static shape, H → H 0 (z; B) as t → ∞, which satisfies whilst MΓ → MΓ 0 as t → ∞ where Γ 0 has a linear profile with slope B/M and mass L over 0 < z < L, For the solution (3.2b) to be valid, it must have Γ 0 ≥ 0 everywhere, or specifically 2M ≥ BL.Indeed, we observe that in simulations with 2M < BL, H and Γ do not approach H 0 and Γ 0 at late times, but rather approach different static solutions: we will discuss this in more detail in § 3.3.Until then, we will focus attention on the case where M is sufficiently large that a solution of (3.2) is approached at late times.
In the surfactant-free problem (Shemilt et al. 2022), the late-time static solution for H also satisfies the ODE (3.2a) but with 2B replaced by B. Hence, the presence of sufficiently strong surfactant effectively doubles the capillary Bingham number in relation to the layer's final shape.In keeping with this observation, figure 3(b) indicates the decreased late-time height of the layer compared with the clean viscoplastic problem. Figure 3(c) shows the evolution of max z (| τw |), the maximum value of the shear stress exerted on the tube wall, showing that it peaks around t ≈ 100 and then approaches B at late times.The comparison with the surfactant-free simulation in the same figure shows that introducing surfactant reduces the peak in the wall stress.

Late-time dynamics of a thin film
To predict the late-time behaviour of the layer, we propose an asymptotic solution for t 1. Inspired by observations from numerical simulations such as figure 3, we make expansions of the form where Γ 0 is given by (3.2b) and H 0 is a solution to (3.2a).For all of the analysis in this section, we will assume that 2M ≥ BL, so (3.2b) is a valid solution.We will also assume that Y −,1 > 0 and Y +,1 < 0 for 0 < z < L, since we observe in numerical simulations (e.g.figure 3) that the whole layer exhibits yielding of type I when it is in this late-time regime.
We will subsequently confirm that the resulting solution is consistent with simulations.
There are two branches of solutions for H 0 , as shown in figure 4(a), which are the same set of static solutions as computed in the surfactant-free problem (Shemilt et al. 2022) but with each corresponding B halved.The upper-branch solutions are strongly deformed, marginally yielded states that the layer approaches at late times.Figure 4(b) illustrates how the deformation and peak height of these solutions is reduced when surfactant is present compared with when the interface is clean, and figure 4(a) quantifies this effect for varying B. The lower-branch solutions in figure 4(a) are unstable near-flat marginally yielded states, which in the clean problem were shown to approximate the minimum initial perturbation required to trigger instability (Shemilt et al. 2022).We will show in § 3.3 that the lower-branch solutions for the surfactant-laden problem have the same significance when M is sufficiently large.Figure 4(a,b) quantify the increased deformation in the lower-branch solutions when surfactant is present compared with when it is not.
Figure 4(c) shows the O(1/t) terms in an example late-time asymptotic solution.There is close agreement between these asymptotics and the numerical solution of the thin-film equations at t = 10 4 .Also, figure 4(c) suggests that Y −,1 ≈ Y +,1 − H 1 , so at late times, the near-wall and near-interface yielded regions have approximately the same size.This implies that the surface velocity, ws , is approximately zero.We investigate this further by plotting the leading-order surface velocity, W 1 , for various M in figure 4(d).As M is increased, MW 1 approaches a fixed curve, indicating that W 1 is approaching zero at a rate proportional to O(1/M).It also shows that W 1 < 0 in all cases, indicating that at late times, surfactant typically induces a reverse flow at the interface.However, at sufficiently large M, the surface velocity is very small and so the interface is essentially immobilised as the layer approaches its late-time configuration.By immobilisation of the interface, we mean that there is no axial surface velocity, but the free surface can still deform and the layer thickness can still evolve.3) for a solution with B = 0.05 and M = 0.5, showing H 1 (solid black), Y −,1 (solid magenta), Y +,1 − H 1 (solid red) and Γ 1 (solid blue).These are compared with the corresponding quantities from a numerical solution of the thin-film equations (2.40)-(2.45) at t = 10 4 , specifically [H(z, t = 10 4 ) − H 0 (z)]Bt (dashed black), Y − (z, t = 10 4 )Bt (dashed magenta), [Y + (z, t = 10 4 ) − H(z, t = 10 4 )]Bt (dashed red) and [Γ (z, t = 10 4 ) − Γ 0 (z)]Bt (dashed cyan) where Γ 0 is defined in (3.2b).(d) Leading-order surface velocity, W 1 , scaled by M, for B = 0.05 and M = {0.25,0.5, 1, 2, 4, 8}.

Effect of varying the Marangoni number on stability and dynamics
To systematically assess the effect of surfactant on the stability and dynamics of the layer, we have run a large set of numerical simulations for various values of B and M, with a fixed initial perturbation to the layer height, A = 0.2. Figure 5 shows the resulting data for the final peak height, max z H(z, t = 10 4 ), which have several characteristics of note.
There is a sharp stability boundary in the data in figure 5.There is a critical value of B, which we call B c , such that for B < B c , there is instability and significant deformation of the layer, whilst for B > B c , there is minimal or no deformation.This stabilisation at high B was identified for the surfactant-free problem (Shemilt et al. 2022), but it can be seen from figure 5 that B c also depends on M. For small M, B c recovers its value for the surfactant-free case but as M is increased, B c decreases towards a different constant value at high M.It is evident that B c also depends on the amplitude of the initial perturbation given to the layer, A, as we discuss in detail in § 3.4 below.Shemilt et al. (2022) showed that B c can be predicted accurately for a large range of A using the marginally yielded static solutions, H 0 (z; B).Here, we find that the same approach can be used to predict B c for large M as well as small M. If we define B m (A) as the value of B such that the solution H 0 to (3.2a) satisfies 1 − H 0 (z = 0; B) = A, then we find B m (0.2) ≈ 0.0289 which can be seen in figure 5 to accurately predict B c at large M. To approximate the stability boundary at small M, we can use the prediction for the surfactant-free problem, which is simply B c ≈ 2B m ≈ 0.0578, since the function H 0 is the same but the corresponding capillary Bingham number for that H 0 is doubled.Figure 5 shows that this predicts the stability boundary well at small M.
To understand why the prediction for large M works, first note that in the surfactant-free problem (Shemilt et al. 2022), the lower-branch static solutions from figure 4(a) correspond to the minimum amplitude of perturbation required to make the layer initially yield, as long as we now also have that MΓ z ≈ B. For large M, surfactant is strong enough that a thin fully yielded region adjacent to the interface develops rapidly before any significant deformation of the layer occurs, so after a rapid adjustment from the initial conditions at early times, we do have MΓ z ≈ B subsequently.Then, whether instability occurs depends only on whether the initial perturbation to the layer height, parametrised by A, was larger than the minimum required to yield, which is represented by the lower-branch static solution for a given B. Or, equivalently, if A is fixed as in figure 5, then instability occurs only if B < B m (A).
In the limit of very strong surfactant, we can derive a simplified evolution equation by exploiting M 1 as a large parameter (see Appendix C).In this large-M asymptotic theory, we find that the surface velocity is zero to leading order, so the interface is immobilised throughout the evolution.Additionally, we find that the motion is governed by the evolution equation, If we replace B with 1 2 B and t with 1 4 t in (3.9), then we exactly recover the evolution equation for the surfactant-free problem ((2.27) of Shemilt et al. 2022).Therefore, solutions for the dynamics with very strong surfactant are the same as solutions for the dynamics without surfactant, but with the capillary Bingham number doubled and time slowed by a factor of four.The four-fold increase in the time scale has been established previously for Newtonian fluids (Otis et al. 1993), but the effective doubling of the yield stress by strong surfactant is (we believe) a new phenomenon.We have already seen this doubling of B in relation to the final shape of the layer in (3.2a) and in the value for B c in figure 5, as described above.However, this theory indicates that for M 1, the doubling of B holds for the entire dynamics, except for a very short period with time scale O(1/M) at the beginning of the evolution when the initially uniform surfactant profile rapidly adjusts to a configuration consistent with the asymptotic theory.
Figure 5 also indicates that, within the unstable region of parameter space, i.e.B < B c , there are two qualitatively different behaviours.At sufficiently large M, there is a region where the final peak height is effectively independent of M.These are the simulations which obey the late-time behaviour described in § 3.2, and so their final shape is the marginally yielded static shape H 0 , a solution to (3.2a), which has no M-dependence.In contrast, for sufficiently small M, not all of the fluid adjacent to the interface fully yields at late times, and so the layer does not enter the late-time regime described in § 3.2 but instead approaches a different final static shape.It can be seen in figure 5 that for small M, the final peak height and, hence, the final shape, depend on both M and B. These results are consistent with the observation that the late-time solution (3.2) exists only for 2M ≥ BL: figure 5 shows there is M-dependence in the final layer shape when 2M < BL.There is, however, also a small band of parameter values in figure 5 where there is still M-dependence in the final shape despite 2M ≥ BL.Although our results suggest that the final shape of the layer in these cases does depend on M, we cannot rule out the possibility that if simulations were run to much longer times, some of these simulations may eventually approach the M-independent final shape, H 0 .Figure 6 shows an example numerical simulation in this small M region.In this simulation, there is a region near z = 0 that exhibits yielding of type II throughout the evolution, even at late times, so there is only partial yielding adjacent to the interface.For this simulation, B = 0.04 and M = 0.08, so 2M < BL.This means that the layer cannot enter the late-time regime described in § 3.2, as then (3.2b) would require Γ to be negative.In figure 5, the boundary between the region with complete surface yielding (simulations that enter the late-time regime from § 3.2) and the small M region with no or partial surface yielding occurs close to the line 2M = BL.Physically, we can interpret the behaviour at small M as the surfactant being too weak for Marangoni effects to be able to fully yield the whole interface.In figure 6, a large portion of the domain does exhibit interface-adjacent yielding, but if M was decreased further, this portion would become smaller and eventually as M → 0, no fully yielded region would exist near the interface, consistent with the behaviour in the surfactant-free case.We note that at late times in the simulation in figure 6, there is complex behaviour near z = 0 with a small region developing where pz changes sign.We find this type of behaviour to be typical in simulations with very small M, but since it occurs at late times when the layer is already near-static, it generally has minimal impact on the global dynamics.Finally, figure 5 demonstrates that surfactant has an appreciable impact on the final configuration of the thin film for M = O(1), corresponding (from (2.38)) to a negligible (O( 2)) reduction in surface tension relative to its mean value.In this thin-film limit, surfactant operates solely through powerful Marangoni effects.and Y + (red), and the thin-film axial velocity, w, represented by the colour map; the middle row shows Γ (green), Γ z (blue) and pz (magenta); and the lower row shows shows the solution in (H pz /B, MΓ z /B)-space, with the black dots corresponding to evenly spaced points along 0 < z < L and the arrows indicating the direction of increasing z.Red diamonds in the first and third panels mark the boundaries of the region where there is yielding at the interface.
3.4.Dependence on the amplitude of initial perturbation When B > 0, the instability is nonlinear and the dynamics of the layer, including its final shape, have some dependence on the size of initial perturbation to the free surface, A. For the surfactant-free problem, Shemilt et al. (2022) showed that the late-time static solutions (solutions to (3.2a) but with 2B replaced by B) can be used both to predict the minimum A required to trigger instability, referred to as A c (B), and to identify a large region of parameter space where the final shape is independent of A. Figure 7(a) illustrates how, in the surfactant-free case, the curve B = 2B m (A) approximates closely the boundary of the region in (B, A)-space where max z H(z, t = 10 4 ) is independent of A. This boundary includes the sharp stability boundary separating simulations where minimal or no growth occurred from the simulations where there is significant growth.As above, B m (A) is defined as the value of B such that 1 − H 0 (z = 0; B) = A, where H 0 are solutions to (3.2a).
Figure 7(b,c) show how introducing surfactant alters this A-dependence, particularly how A c is increased for larger values of M. For M = 0.6 (figure 7b), A c is increased for each value of B compared with the surfactant-free case.When surfactant strength is increased again to M = 6 (figure 7c), A c is further increased and, analogously to the surfactant-free case, the curve B = B m (A) predicts the entire boundary of the region where the final shape is independent of A, including predicting A c .The lower branch of static solutions in figure 4(a) correspond to the values of B m that predict A c in figure 7(c), illustrating that these lower-branch solutions correspond to the minimal initial perturbation to the free surface required to trigger instability when surfactant is strong.The accuracy of B = B m (A) in predicting the whole boundary of the A-independent region in figure 7(c) provides further evidence of the effective doubling of B by introducing strong surfactant compared with the surfactant-free case, and highlights that the static solutions, H 0 (z; B), can provide significant insight into the stability and dynamics of a layer with strong surfactant.Figure 7 also illustrates that, while there is non-trivial dependence of A c on M, accurate upper and lower bounds for A c for any M are provided by the curves B = B m (A) and B = 2B m (A), respectively.

Thick films and liquid plug formation
Whilst the relative simplicity of the thin-film equations allows for more detailed analysis, liquid plug formation cannot occur in that theory.To assess the effect of surfactant on the dynamics leading to plug formation, we now consider the long-wave equations (2.27)-(2.33),which model the evolution of layers which are not thin.The layer thickness is described by the parameter , which is the ratio of the average thickness to the tube radius after a small adjustment to account for the change in volume induced by having a finite amplitude initial perturbation (see (2.48)).
The minimum volume of fluid required to form a liquid plug corresponds to a layer thickness of ≈ 0.107 for the length of domain we are using (Everett & Haynes 1972).Gauglitz & Radke (1988) conducted the first numerical simulations of the Newtonian problem without surfactant, which predicted a slightly larger critical thickness, crit ≈ 0.12, due in part to the restriction of only being able to run simulations to a relatively short, finite time.Viscoplastic rheology can significantly increase crit when the capillary Bingham number is sufficiently high (Shemilt et al. 2022).Here, we investigate how crit is affected by the additional presence of surfactant.
We run all simulations to time t = 10 4 , which is an order of magnitude longer than in the investigation by Shemilt et al. (2022), to capture the dynamics around plug formation as accurately as possible.Following the approach of many previous studies, we stop the simulations early if min z R ≤ 0.3, or equivalently if max z H ≥ 0.7/ , as this has been found to indicate that a plug is imminently about to form (Halpern & Grotberg 1992, 1993;Halpern et al. 2010;Shemilt et al. 2022).We then identify the time at which the simulation is stopped as the plugging time, t p .We will describe the thick-film results in this section using the unscaled capillary Bingham and Marangoni numbers, B and M, defined in the long-wave theory (2.38), rather than the scaled versions that arise in the thin-film theory.
One can convert back to the scaled, thin-film parameters by dividing B and M by 2 .
4.1.Example evolution showing plug formation Figure 8(a) shows a numerical solution of the long-wave equations.The layer thickness is = 0.14, meaning that there is sufficient volume of fluid that a plug could form, and it can be seen in figure 8(b) that one does form at around t ≈ 268.Around 267 t 268, the layer grows rapidly towards the centre of the tube, indicating that it is about to form a plug when the simulation is stopped.Introducing surfactant to a Newtonian layer with comparable thickness has been shown to increase the plugging time approximately four-fold (Ogrosky 2021).We find that adding surfactant to a viscoplastic layer can increase the plugging time by a much greater amount.In the example shown in figure 8(b), plug formation in the surfactant-laden viscoplastic film takes approximately eight times longer than when there is no surfactant, which in itself takes longer than plugging in the surfactant-free Newtonian simulation.
Figure 8(a) shows that by t = 190, the layer is exhibiting yielding of type I everywhere, similar to the typical thin-film behaviour seen in figure 3(a).However, at later times, a region with type IV yielding, where the fluid is yielded near the interface but rigid adjacent to the wall, appears near z = 0 and then expands towards the centre of the domain.The flow in this region is in the negative z-direction, opposing the flow in the rest of the layer.However, the velocity in this region is small so we expect any effect on the global dynamics to be minimal.This behaviour is similar to what is observed in the long-wave problem without surfactant, but in that case, the region near z = 0 is fully rigid (Shemilt et al. 2022).Whilst we expect this long-wave theory to be a good approximation to the overall dynamics of the thick film, small details such as this behaviour near z = 0 would need to be confirmed with full two-dimensional simulations.The more profound effect of surfactant that can be seen in figure 8(a) is the development of a sizeable fully yielded region near the interface all along the layer, where the shear is opposing the general flow and so delaying the formation of a plug.

Thick-film dynamics at large Marangoni numbers
To explore the effect of increasing the surfactant strength on the evolution of a thick film, we propose an expansion for M 1, Since the surface stress must be finite, Conservation of the total amount of surfactant in the long-wave theory then implies The surfactant transport equation (2.30) at leading order in M, after rearranging and using (4.2), gives Therefore, unlike for the thin-film case (Appendix C), the surface velocity is not zero for large M, in general.From (4.3), the surface velocity can be deduced from the difference between the local time derivative of the film thickness, ∂ t R 0 , and the globally averaged rate of change of the thickness, via the integral of ∂ t R 0 .Non-zero W 0 leads to small gradients in surfactant concentration.The results (4.2) and (4.3) can be understood physically as follows: the surfactant is so strong that any local change in the shape of the interface induces a surface velocity that redistributes surfactant so that the concentration remains (almost) globally uniform, and if the total surface area of the interface has changed, this will also induce a change in the mean concentration.We have appealed neither to the equation of state for surface tension (2.19), except that σ (Γ ) / = 0, nor to the liquid's rheology to reach (4.2) and (4.3).Hence, within long-wave theory, the results (4.2) and (4.3) are generic for any type of fluid and any equation of state that has σ (Γ ) / = 0. Figure 9 compares a solution of the full long-wave system (2.27)-(2.33)with M = 10 with the results (4.2) and (4.3).It shows good agreement between Γ and w s from the numerical solutions and the leading-order approximations (4.2) and (4.3), despite the value of M not being extremely large.The spatial gradients in Γ (figure 9b) are small in the numerical solution, and we expect these would become smaller if the surfactant strength was further increased, in line with the large-M theory.Figure 9(c) shows that the surface velocity induced as plug formation occurs is in the negative z-direction, and is strongest immediately before the simulation is stopped.This is because there is a rapid decrease in the local surface area of the interface near z = L during the fast growth before plug formation, which induces a surface velocity to redistribute surfactant across the layer.These are compared with G 0 (t) (dashed), the approximation from the large-M asymptotic theory, which is evaluated via (4.2) using the numerical solutions from panel (a) as proxies for R 0 .(c) Surface velocity, w s (solid), at the same time points.These are also compared with the corresponding approximation from the large-M theory (4.3).
The results from figure 9 suggest that the large-M theory can be used to better understand the behaviour of thick films when surfactant is strong.However, we have assumed a simple linear equation of state for surface tension (2.19), which limits how large we are able to make M in simulations while retaining accuracy.Figure 9(b) suggests that when the simulation is stopped, the increase in Γ from its original value was approximately 3 %.We have found that this is approximately independent of M, so for values of M close to or larger than M ≈ 30, we expect the theory to break down since σ will approach zero as a plug forms.For this reason, M = 10 is the largest Marangoni that we have used in the long-wave simulations.A more complex nonlinear equation of state would have to be used to access higher values accurately.

Critical layer thickness for plug formation
Figure 10(a) compares the computed critical layer thickness for plug formation to occur, crit , when surfactant is present with its value when the interface is clean.Both with and without surfactant, increasing B induces an increase in crit .The increase is small when B is small, but is significant at larger B. When surfactant is present, the increase in crit is initiated at B ≈ 10 −3 , compared with at B ≈ 2 × 10 −3 when the interface is clean.Beyond these values of B, the increase is amplified by the presence of surfactant, such that the value of crit is approximately 30 % larger than when the interface is clean.Figure 10(b) illustrates the effect on the dynamics of varying M and .It shows data from many simulations, some of which are stable and some form plugs.The boundary between the stable and unstable regions in the data can be identified as crit , so the dependence of the critical layer thickness on M can be seen.At low M, the critical thickness for the surfactant-free problem is recovered, and as M is increased, crit increases to a significantly higher value for large M. Therefore, for this value of B, there is a range of layer thicknesses 0.14 0.185 where increasing the surfactant strength sufficiently, without changing any other parameters, can suppress plug formation.This stabilisation of the system by increasing the surfactant strength resembles what was seen in the thin-film system (figure 5).Again, surfactant has a powerful effect at low concentrations: the range 10 −2 < M < 10 −1 , over which crit changes appreciably in figure 10 1 % and 10 %. Figure 10(b) also quantifies how the plugging time increases as M is increased.We know increasing B extends the plugging time (Shemilt et al. 2022), and figure 10(b) shows it can be extended further again by introducing surfactant.However, the profound stabilisation effect by surfactant is not due mainly to slowing of the dynamics as it is in the Newtonian problem (Halpern & Grotberg 1993).Instead, figure 10(b) highlights how the increase in crit due to rigidification and stabilisation of the layer by yield stress is amplified by the presence of sufficiently strong surfactant.
The value of crit and the dynamics leading to plug formation are necessarily dependent on the initial perturbation applied to the layer, which here is parametrised by A. With a larger initial perturbation, there can be more yielding at larger values of B. However, we have found that varying the initial conditions does not qualitatively affect the results presented.The finite time to which we run simulations, t = 10 4 , will have some effect on the computed value of crit , but figure 10(b) shows that the plugging time increases rapidly close to the computed value of crit , suggesting that there is good convergence to the true value of crit .

Discussion
We have developed two models, using long-wave and thin-film theories, of the capillary instability of a viscoplastic layer coating a cylindrical tube where insoluble surfactant is present at the air-liquid interface.This flow can be strongly modified both by the strength of the yield stress, described by the capillary Bingham number, and the strength of the surfactant, described by the Marangoni number.We showed that the layer can exhibit five qualitatively different types of yielding, which are different combinations of fully yielded shear flow, pseudo-plugs and rigid plugs, depending on the relative strengths of the capillary and Marangoni forces.Numerical simulations highlighted the complex dynamical transitions between these yielding types that can occur as the layer evolves.
For thin layers, we quantified how introducing surfactant enhances the stabilising effect of the yield stress and induces an effective doubling of B when M is sufficiently large.Asymptotic analysis of the late-time behaviour in § 3.2, valid for moderate and large M, showed that the final static, marginally yielded configuration of a surfactant-laden layer coincides with that of a surfactant-free layer but with B exactly doubled.Static, marginally yielded solutions were also shown to accurately predict the critical capillary Bingham number above which instability is suppressed, for both small and large Marangoni numbers (figure 5).When M is asymptotically large, the entire thin-film dynamics coincides with that of a surfactant-free layer but with time slowed by a factor of four and B doubled.
For thicker layers, we quantified how the known effect of yield stress in delaying or suppressing plug formation (Shemilt et al. 2022) is amplified by introducing surfactant.Using numerical solutions of the long-wave equations, which describe the motion of thick films, we showed that the approximately four-fold increase in plugging time when strong surfactant is introduced to a thick Newtonian layer (Otis et al. 1990;Halpern & Grotberg 1993) can become much larger when the liquid is viscoplastic (figure 8).The critical layer thickness required for a plug to form is also known to increase as the capillary Bingham number is increased (Shemilt et al. 2022).In figure 10, we quantified how this increase is amplified when sufficiently strong surfactant is present; the threshold for surfactant to have an appreciable effect corresponds to a relative surface-tension reduction of only a few percent.By examining the limit of large Marangoni number in § 4.2, we found that surface velocities are induced by local changes to the shape of the interface during plug formation, which act to redistribute surfactant and uniformly raise the concentration globally.
The results suggest a mechanism for the stabilising effect of pulmonary surfactant in the small airways, particularly in diseases where mucus yield stress is increased such as cystic fibrosis (Patarin et al. 2020).Figure 8(b) illustrates the delay to plug formation caused by surfactant.In these simulations, plug formation occurs at t p ≈ 35 when there is no surfactant and t p ≈ 268 when there is surfactant.We can relate these dimensionless times to real time scales in the lungs by taking the following physiologically relevant parameters values for a 12th generation airway: airway radius, a = 0.4 mm (Hsia, Hyde & Weibel 2016); equilibrium surface tension, σ 0 = 30 mN m −1 (Chen et al. 2019) and viscosity, η = 0.01 Pa s (Lai et al. 2009).Using these values, t p ≈ 35 and t p ≈ 268 correspond to real plugging times of t * p ≈ 1.7 s and t * p ≈ 13 s.Whilst there is significant variation in measurements of rheological parameters for mucus, this suggests that for a 12th generation airway, the plugging time in the surfactant-free simulation is less than the time of a breathing cycle and it is longer than a breathing cycle in the surfactant-laden simulation.This provides evidence that surfactant can make plug formation less likely to occur in an airway.
The results presented in figure 10 for the critical layer thickness for plugging to occur also have physiological relevance.Measurements of healthy pulmonary surfactant suggest that Marangoni numbers around M ≈ 1 are likely to be observed (Schürch et al. 2001), although there is significant variation in measured values.Figure 10(b) illustrates how the critical thickness is increased when M ≈ 1 compared with when M is low.There is clinical evidence that in lungs affected by cystic fibrosis, surfactant function is impaired (Gunasekara et al. 2017) and that this is linked to increased prevalence of airway obstructions (Griese et al. 2004).Our results are consistent with these clinical observations, and suggest a mechanism for how surfactant deficiency may destabilise airways.Moreover, as discussed by Shemilt et al. (2022), mucolytic therapies commonly used in cystic fibrosis can significantly decrease mucus yield stress (Patarin et al. 2020), which could destabilise airways and trigger plug formation if the mucus layer is thick.One of these mucolytic drugs, rhDNAse, has been linked with increased lung exacerbations and a decline in lung function in patients with idiopathic bronchiectasis (O'Donnell et al. ).The evidence of surfactant amplifying yield-stress effects suggests surfactant needs to be accounted for when considering the impacts of mucolytic drugs that significantly lower the mucus yield stress.
Our results also suggest that introducing surfactant can decrease the peak value of the shear stress induced on the tube wall during the evolution of the liquid layer (figure 3c).This suggests that surfactants in airways are likely to provide protection against epithelial cell damage, which can occur when a sufficiently large stress is exerted on the airway wall (Huh et al. 2007).This is consistent with previously reported effects of introducing surfactant to a Newtonian layer (Romanò et al. 2022).
The theory we have used here could be easily modified to study other related flows.For example, the surface Marangoni force could be replaced by a force caused by an oscillating air flow, something that has been studied previously in the case that the liquid is Newtonian (Halpern & Grotberg 2003).Beyond the application to airway modelling, the present theory could be used to study the effect of surfactant on other capillary flows of viscoplastic fluids.For example, one possible industrial application that has received recent modelling attention is ink-jet printing (Jalaal et al. 2021;van der Kolk et al. 2023), where the effects of surfactant have not yet been studied but may be important.The flow region maps that we have introduced (figures 2, 3a and 6) may also prove useful in other problems by illustrating the transitions between different types of yielding.
There are many limitations to the modelling approach we have taken.Several physiologically relevant effects have not been included in the model, such as gravity, air flow and more complex liquid rheologies.However, the relative simplicity of the model has allowed for extensive exploration of parameter space to systematically examine the effect of surfactant on this flow.We have assumed a linear relation between surface tension and surfactant concentration, despite this typically being nonlinear for pulmonary surfactants (Schürch et al. 2001).With this assumption, the model is simplified but still retains sufficient complexity to capture the key effects of surfactant on the flow.As discussed in § 4.2, whilst we have explored a wide range of M, the linear surfactant equation of state limits how large we can make M in our long-wave simulations.To fully explore the large-M limit in the future, a nonlinear equation of state should be used.The quasi-one-dimensional long-wave theory has been shown to provide a good approximation to the dynamics of thick films in similar problems (e.g.Halpern et al. 2010) so we expect it to perform well here also, but our predictions await validation against fully two-dimensional CFD simulations.Long-wave theory is known to break down immediately before plug formation (Johnson et al. 1991), so to extend the analysis to study coalescence and the post-coalescence phase would require a fully two-dimensional model, in line with what has been done in the Newtonian case (Romanò et al. 2022).
In this study, we have quantified how surfactant can amplify the effect of viscoplastic rheology during the capillary instability of a liquid film coating a tube.In addition to slowing the dynamics, sufficiently strong Marangoni forces can induce an effective doubling of the capillary Bingham number for thin films.For thick films with a large enough capillary Bingham number, the critical thickness required for plug formation to occur can be significantly increased by the presence of surfactant.These results suggest a mechanism for how pulmonary surfactant can stabilise airways and how surfactant deficiency can contribute to the prevalence of mucus plugging in obstructive lung diseases.

971Figure 1 .
Figure 1.(a) Sketch of the model geometry.The air-liquid interface is located at r = R(z, t).Insoluble surfactant is present at the interface, with (non-dimensionalised) concentration Γ .(b) Illustration of a possible axial velocity profile in the liquid layer, w.Fully yielded, shear-dominated regions (white) are shown adjacent to the interface (R ≤ r ≤ Ψ − ) and adjacent to the wall (Ψ + ≤ r ≤ 1), with a plug-like region (grey) in between (Ψ − < r < Ψ + ).
Figure 2(c) illustrates where in (H pz /B, MΓ z /B)-space each of types I-V occurs.The picture is similar to figure 2(b), which was described in detail above, but is somewhat simplified.
(a)  for 80 ≤ t ≤ 130.Once the shock waves have propagated to the side boundaries, the evolution enters a late-time regime where the upper and lower yielded regions extend across the whole 971 A24-14 https://doi.org/10.1017/jfm.2023.

Figure 3 .
Figure 3. (a) Snapshots from a numerical solution of the thin-film evolution equations (2.40)-(2.45)with B = 0.04, M = 0.2, A = 0.2 at t = {0, 5, 80, 100, 130, 250, 500, 2000}.At each t, there are three panels: the top panel shows the layer evolving, with Y − (cyan) and Y + (red), and the thin-film axial velocity w represented by the colour map; the middle panel shows plots of pz (magenta), Γ (green) and 10Γ z (blue); the bottom panel shows the solution in (H pz /B,MΓ z /B)-space, in dotted black lines, with the dots corresponding to points evenly spaced along the domain 0 < z < L and the arrows indicating the direction of increasing z.Red diamonds in the first and third panels mark the boundaries of the region where there is yielding at the interface.(b) Time evolution of max z H from the same simulation (solid) compared with the evolution of max z H from a Newtonian surfactant-laden simulation with (B, M) = (0, 0.2) (dashed), a surfactant-free viscoplastic simulation with (B, M) = (0.04, 0) (dash-dotted) and a surfactant-free Newtonian simulation with (B, M) = (0, 0) (dotted).(c) Time evolution of max z Y − (solid red), max z (H − Y + ) (solid blue) and max z (| τw | − B) (solid black) for the simulation in panel (a).Also shown are plots of max z Y − (dash-dotted red) and max z (| τw | − B) (dash-dotted black) for the surfactant-free viscoplastic simulation (B = 0.04, M = 0).

Figure 5 .
Figure 5. Data from thin-film simulations with A = 0.2 and various B and M. Each coloured point corresponds to one simulation with the colour indicating the final peak height.Contours interpolated from the same data are also plotted in black, which are evenly spaced and in the range 2.4 ≤ max z H ≤ 2.8.The red lines are 2M = BL (solid), B = B m (A = 0.2) ≈ 0.0289 (dashed) and B = 2B m (A = 0.2) ≈ 0.0578 (dash-dotted).

971Figure 6 .
Figure6.Snapshots from a thin-film simulation with B = 0.04, M = 0.08, A = 0.2, at t = {70, 120, 1100, 9000}.The upper row of panels shows the layer height evolving, with Y − (cyan) and Y + (red), and the thin-film axial velocity, w, represented by the colour map; the middle row shows Γ (green), Γ z (blue) and pz (magenta); and the lower row shows shows the solution in (H pz /B, MΓ z /B)-space, with the black dots corresponding to evenly spaced points along 0 < z < L and the arrows indicating the direction of increasing z.Red diamonds in the first and third panels mark the boundaries of the region where there is yielding at the interface.

971Figure 7 .
Figure 7. Data from thin-film numerical simulations with (a) M = 0, (b) M = 0.6 and (c) M = 6.Each coloured dot corresponds to one simulation, with the given values of A and B, where the colour indicates the final maximum height of the layer.The same data for max z H(z, t = 10 4 ) are linearly interpolated and plotted as black contour lines.The two magenta curves in each panel are B = B m (A) (solid) and B = 2B m (A) (dash-dotted).

971Figure 8 .
Figure 8.(a) Numerical solution of the long-wave equations (2.25)-(2.33)with A = 0.2, B = 0.001, M = 0.02 and = 0.14, showing the transition towards plug formation.The upper row shows the evolution of the layer height, with Y − (cyan) and Y + (red) also shown, and the colour corresponding to the magnitude of the axial velocity, |w|.The lower row shows p z (magenta), Γ z (blue) and Γ (green).(b) Time evolution of max z H for the same simulation, compared with max z H from simulations with (B, M) = {(0, 0), (0.001, 0), (0, 0.02), (0.001, 0.02)}, showing the combined delay to plug formation by surfactant and yield stress.

971Figure 9 .
Figure 9. Results from a numerical solution of the long-wave equations with = 0.14, M = 10, B = 0.001 and A = 0.25.(a) The interface position, R, shown at various time points.The last time point is t = t p = 410.69,when the simulation is stopped.(b) Surfactant concentration, Γ (solid), at the same time points as in panel (a).These are compared with G 0 (t) (dashed), the approximation from the large-M asymptotic theory, which is evaluated via (4.2) using the numerical solutions from panel (a) as proxies for R 0 .(c) Surface velocity, w s (solid), at the same time points.These are also compared with the corresponding approximation from the large-M theory (4.3).

Figure 10
Figure10.(a) Critical layer thickness, crit , as a function of B, for a surfactant-free layer (M = 0) and a layer with surfactant (M = 0.4).All simulations have A = 0.25.The value of crit computed is such that a simulation with = crit + 0.001 forms a plug before t = 10 4 and a simulation with = crit − 0.001 has not formed a plug by t = 10 4 . (b) Data from long-wave simulations at various values of and M, with A = 0.25 and B = 0.0024.Grey crosses indicate simulations where a liquid plug formed, while black dots indicate simulations where a plug did not form.Within the plugging region, the grey contours indicate the plugging time, t p .