Trapping and displacement of liquid collars and plugs in rough-walled tubes

A liquid film wetting the interior of a long circular cylinder redistributes under the action of surface tension to form annular collars or occlusive plugs. These equilibrium structures are invariant under axial translation within a perfectly smooth uniform tube and therefore can be displaced axially by very weak external forcing. We consider how this degeneracy is disrupted when the tube wall is rough, and determine threshold conditions under which collars or plugs resist displacement under forcing. Wall roughness is modelled as a non-axisymmetric Gaussian random field of prescribed correlation length and small variance, mimicking some of the geometric irregularities inherent in applications such as lung airways. The thin film coating this surface is modelled using lubrication theory. When the roughness is weak, we show how the locations of equilibrium collars and plugs can be identified in terms of the azimuthally averaged tube radius; we derive conditions specifying equilibrium collar locations under an externally imposed shear flow, and plug locations under an imposed pressure gradient. We use these results to determine the probability of external forcing being sufficient to displace a collar or plug from a rough-walled tube, when the tube roughness is defined only in statistical terms.


I. INTRODUCTION
Surface-tension driven flows in liquid-lined lung airways are of physiological importance in driving airway closure (through the Rayleigh instability, causing redistribution of airway liquid into occlusive liquid plugs, and through compressive loading of capillary forces on the flexible airway wall) and subsequent airway opening (by displacement of liquid plugs and airway inflation). Significant insight into the physical processes underlying mechanisms of airway closure and reopening have come from studies of idealised model problems, many of which have wider relevance to two-phase flow in porous media and microfluidics. In particular, models based on lubrication theory of the initial Rayleigh instability in a uniform liquid-lined tube, describing the formation of annular collars [17] and liquid bridges (or plugs) [13], have been extended to account for numerous features of the complex in vivo airway environment in health and disease, including factors such as wall elasticity and airway collapse, surfactants, imposed shear due to airflow, gravity, airway centreline curvature and non-Newtonian rheology (reviewed in [2, 14-16, 19, 20, 27]).
The present study addresses a curious aspect of the original problem studied by Hammond [17]. A thin, initially uniform liquid film coating the interior of a long circularly-cylindrical tube, subject to no-flux conditions applied at either end of the tube, redistributes under the action of surface tension into a set of annular collars, connected to each other by a slowly draining film that remains continuous in the absence of evaporation and destabilising disjoining pressure. While the collar shape is well approximated as a surface of uniform curvature (an unduloid that meets the tube wall with zero contact angle [7]), the collar location is not so readily determined. In practice, collars will either stabilize at either end of the tube (with centres coincident with no-flux boundaries), or they may migrate towards a boundary, even reversing direction over very long time-scales under the influence of very small differences in the draining flows near each edge of the collar [29]. This sensitivity reflects a degeneracy in the underlying model, whereby an equilibrium wetting collar (or, if sufficient fluid is available, an occlusive liquid plug) can in principle be located anywhere along a sufficiently long uniform circular cylinder.
In reality, lung airways (and other tubes arising in natural environments) are not perfectly uniform cylinders. Major geometrical imperfections (such as centreline curvature) have already been shown to disrupt collars sufficiently to form axial rivulets [18,25]. Here we address a more subtle distortion, imagining that the tube has small random perturbations to its shape (in the case of an airway, arising for example from protruding epithelial cells, inhaled debris or mucosal buckling), which are described by a Gaussian random field of prescribed variance and correlation length. We seek the conditions under which such perturbations are sufficient to stabilize collars at discrete locations along the tube. Small shape perturbations represent a singular limit of Hammond's problem, regularizing the degeneracy associated with axial collar translation.
We characterise the wall shape in statistical terms (representing the uncertainties and natural variability that are inherent to physiological systems) and correspondingly express outcomes in probabilistic terms. We build on prior deterministic studies of the influence of axially periodic axisymmetric corrugations along the tube on collar formation, addressing linear [46] and weakly nonlinear [47] stability as well as the fully nonlinear dynamics leading to tube occlusion [1,44], which typically show alignment of near-equilibrium structures with tube constrictions. Non-axisymmetric perturbations to the exterior of a liquid-lined cylinder can also influence film distributions [28]. Our study complements the extensive literature on liquid plugs in lung airways, which focuses primarily on plug displacement under forcing, with studies focusing on the role of inertia [9], wall flexibility [22,50], surfactant [10,45], dynamical effects including plug rupture [12,23,32,42], gravity [40,51], interactions with bifurcations [43,52] and yielding behaviour [24,49].
This study extends this prior work by addressing the role of random geometric imperfections on collar and plug dynamics. For tubes having weak non-axisymmetric roughness, we show how equilibrium collar and plug locations are defined in terms of the azimuthallyaveraged wall shape, motivating the study of axisymmetric tubes with axially non-uniform shapes. For such tubes, having randomly distributed rather than periodic constrictions, we derive algebraic conditions for the existence of stable capillary equilibria. We use two interpretations of stability: one involving perturbations driven by surface tension effects alone; the other involving an imposed external perturbation (here we consider shear from a core flow displacing collars, or an external pressure gradient displacing occlusive liquid plugs). For applications in which the tube shape is defined only at a statistical level, we then determine the probability that, over many realisations, a given external forcing is sufficient to displace an isolated collar or plug. We first present results for collars (Secs. II-V) and briefly discuss extensions to liquid plugs in Sec. VI below, giving displacement probabilities in an explicit analytic form.

II. MODEL FOR SLENDER LIQUID COLLARS
We consider a rigid hollow cylinder with a rough interior wall coated by a thin layer of fully-wetting Newtonian liquid of viscosity µ and surface tension σ ( Figure 1). We assume the wall roughness amplitude and the liquid-layer thickness are of comparable magnitude and that both are much smaller than the cylinder's mean internal radius a 0 . We express the liquid-layer thickness as a 0 h(θ, z, t) and the internal radius of the cylinder as a 0 (1− ηa(θ, z)), where a is a realisation of a Gaussian random field having zero mean and O(1) exponential covariance. Here 1 is the ratio of the spatially-averaged liquid-layer thickness to a 0 , η is the ratio of the roughness amplitude relative to the liquid-layer thickness, θ ∈ [0, 2π) measures the azimuthal angle around the tube, a 0 z measures axial distance and time t is scaled on a 0 µ/ 3 σ. The tube has length a 0 L. The liquid layer is subject to an imposed axial shear stress 2 στ /a 0 (due to flow of gas in the core of the tube) and satisfies no-slip and no-penetration conditions at the cylinder wall. We neglect gravity and inertia and assume the cylinder's length a 0 L significantly exceeds its radius.
In the leading-order lubrication approximation, the dimensionless liquid-layer thickness h(θ, z, t) satisfies [25,29] where subscripts denote derivatives. The film adjusts under the action of the imposed shear stress τ and pressure gradients arising from the non-uniform curvature of the gas-liquid interface ( Figure 1). The function g(θ, z, t) measures the distance of this interface from the mean wall location; in the expression for the linearized mean curvature in (II.1b), the term g represents the azimuthal curvature that drives collar formation. We impose periodic . .
z (Axial periodicity is imposed for computational convenience but later we will relax this condition.) We impose h = 1 at t = 0, so that V 0 = 2πL. Since disjoining pressure effects are neglected, the film remains continuous with h > 0 everywhere, although it can become very thin over much of the domain at large times.
In order to be compatible with the periodic boundary conditions, we consider the wall roughness field a(θ, z; ω) to be a doubly periodic stationary Gaussian random field with zero mean and covariance function where r 1 and r 2 are O(1) dimensionless correlation lengths associated with the θ and z directions, comparable in magnitude to the cylinder radius (min(r 1 , r 2 ) η ensures that the assumptions of lubrication theory are not violated). The sin 2 form in each direction arises from sampling a planar 2D field with exponential covariance around a circle [35]. Here we allow r 1 and r 2 to be varied independently; our approach below readily accommodates other choices of k, such as which represents the restriction of a 3D isotropic Gaussian field with correlation length r to a long cylinder of unit radius. ω identifies a(θ, z; ω) as a sample drawn from a probability distribution and it implicitly labels realisations of the model, although we largely suppress it in what follows. We use a Karhunen-Loéve decomposition to generate samples of a numerically [30]. The associated film distribution is obtained by solving (II.1) using a finite-element method, implemented in COMSOL Multiphysics, with care taken to resolve fine-scale structures that emerge at large times. A realisation of a(θ, z), and its azimuthal average a(z), are shown in Figure 2(a,b), where the azimuthal averaging operator is defined as The interfacial area (equivalent to a capillary surface energy) associated with (II.1) can be expressed (to the appropriate level of accuracy) as Imposing periodic boundary conditions, integration by parts demonstrates that Imposed shear injects energy into the system, such that A can be expected to evolve to a statistically steady state at large times. However in the absence of shear, (II.7) shows that A t ≤ 0, so that A continually diminishes as the film evolves, although typically the system approaches a metastable state that does not represent the global energy minimum of the system. In a perfectly cylindrical tube (η = 0), the global minimum has the entire fluid volume confined within a single equilibrium collar of the form where the collar may lie anywhere in the domain; elsewhere h is vanishingly small but nonzero. The collar described by (II.8) has zero contact angle at each of its effective contact lines (h 0 = h 0z = 0 at z = z ± π ≡ z ± ). The corresponding energy A is independent of z to this level of approximation; in practice the collar location becomes sensitive to fine details of the external film distribution [29].
The simulation in Figure 2 illustrates the redistribution of a film within a rough-walled tube in the absence of external shear. In this example the initially uniform film evolves into two collars separated by regions in which the film becomes very thin. The wall roughness is isotropic (Figure 2a), with the azimuthally averaged roughness having two prominent maxima within the domain (Figure 2b) on which collars form (Figure 2c). The film remains continuous, with fluid continually draining from the ultra-thin film into the neighbouring collars. The pressure is low and uniform within each collar (Figure 2d), rising abruptly across the effective contact lines, which are transverse to the tube axis and spaced a distance 2π apart across each collar. Unlike the situation in a perfectly smooth tube, for which collars can migrate axially over long time intervals [29], here collars are pinned by the underlying topography. Despite significant azimuthal variation in wall roughness, surface tension drives the liquid towards a more axisymmetric film distribution. The partitioning of the initial fluid volume between the two collars depends on the initial conditions.
As we illustrate below, sufficiently large shear displaces the collars, ultimately leading to travelling-wave solutions for which collars migrate repeatedly through the domain (because of the periodic boundary conditions). We identify the loss of stability of all stationary collars to travelling-wave states as the signature of the condition under which liquid can be permanently displaced from a tube, under suitable outlet conditions.
Below we seek to identify the conditions defining the location of collars (such as those in Figure 2) and their stability to imposed shear. In particular, we wish to describe τ c (ω), the largest shear for which a particular tube realisation a(θ, z; ω) can support at least one stationary stable collar. From this we wish to deduce the collar displacement probability P(τ ). For a tube drawn randomly from a sample with specified covariance (II.3), this is the probability that τ c is below a given shear τ , This quantity indicates the likelihood that, over multiple realisations, τ is sufficiently large to displace all stationary collars of volume V 0 or less that may form in a tube with geometric features characterised by η, r 1 , r 2 and L. A sufficiently large value of τ , ensuring a value of P approaching unity, can be used as a criterion to ensure reliable removal of liquid from a rough tube under imposed shear.
We first address the simplified case when the roughness is axisymmetric (Sec. III), and then return to the non-axisymmetric case in Sec. IV. In Sec. V we determine τ c (ω) by direct simulation and P(τ ) by a Monte-Carlo method, supplementing these predictions with asymptotic approximations assuming weak roughness (η 1). These yield explicit predictions showing how P(τ ) is related to the governing parameters, which we extend to liquid plugs in Sec. VI.

A. Collar locations under zero shear
For τ = 0, at large times there is a weak flux through the effective contact lines of the collar, driven by slow draining from the adjacent film. We address here the equilibrium structure itself, enforcing uniform pressure in the collar and imposing zero effective contact angle at each contact line. Thus we solve The collar solution may be written where A, z , z ± and P are to be determined. The four contact-line conditions plus the volume constraint It is helpful to simplify this problem by assuming 0 < η 1. For η = 0, z is undeter-

8). A regular expansion of
A, P and z ± in powers of η, using gives from the contact-line conditions a ± + A 1 + P 1 = 0 and (V /4π 2 )z ± 1 = a ± z , where a ± ≡ a(z ± π). The volume constraint implies P 1 = −(1/2π) z +π z −π a dz. Thus z is determined by the condition that the contact lines meet the wall at identical radial distances from the tube centreline, For a rough wall of length L this condition identifies a finite number of possible collar locations, a subset of which are realised in practice. The collar width becomes implying for example that a collar that sits on a constriction, for which a − z > 0, a + z < 0, is shortened, whereas a collar that sits in a dilation is elongated. An energetic argument involving the interfacial energy (II.6) shows how stable collars cannot be greater than 2π in length (Appendix A). This is supported by numerical evidence (such as Figure 3a) that collars sitting on constrictions are stable under zero shear; we found no evidence of stable collars sitting in dilations under zero shear.
The stability of a collar under zero shear relies on the interfacial (capillary) energy (II. 6) being at a local minimum. This is distinct from the stability of such collars to weak forcing by shear, which we now consider.

B. Stationary collars under shear
To interpret numerical observations in Figure 3 for τ > 0 we can again use a perturbation analysis assuming both η and τ are small, using the single-collar solution (II.8) as a starting point. As suggested by Figure 3(b), we anticipate that the collar is connected to an external thin film through thin transition regions at each contact line. Within each transition region, we expect the dominant balance of terms in (III.1) to be where q is the local volume flux. Solutions of this inner problem must match the collar (for We can then seek a uniformly steady solution of (III.1).
We therefore develop an expansion in τ 1, formally assuming τ and η are of comparable magnitude as τ → 0. Details of asymptotic matching are provided in Appendix B; here we summarise the conditions determining the collar shape and location. Since the ODE in each transition region (III.5) gives rise to logarithmic terms in the far-field, we describe the collar in |z − z | < π using the expansion Here h 0 is given by (II.8) and z denotes the collar mid-point at which Within the collar, the dominant balance of terms in (III.1) becomes The flux through the collar is neglected in this approximation but shear induces a pressure gradient associated with an internal recirculating flow, visible in Figure 3 The leading-order shape satisfies p 0z = 0, so that h 0 is given by At the following order, p 2z = 3/(2h 0 ), describing the stirring flow in the collar generated by external shear. Integrating, where A 2 and B 2 are constants of integration and F is an odd function satisfying F + F = We use (III.12) to consider the stability of collars in the presence of shear. For a given wall shape over a domain of length L, we anticipate a set of discrete solutions, giving z 0 in terms of the parameter Ξ ≡ τ /(V η). Displacement of the collar due to changes in this parameter are given to leading order by For a − z − a + z > 0 (< 0), as is the case for a collar sitting on a constriction (in a dilation), the collar is displaced downstream (upstream) by increasing shear, which we conjecture implies stability (instability) under shear. This conjecture is consistent with observations in Figure 3.
The threshold condition at which a collar loses stability, which corresponds to a saddle-node bifurcation when collar locations are mapped against Ξ (as illustrated below), is therefore (III.12) subject to a − z = a + z . In other words, the limiting collar locations are defined by points on the wall, an axial distance 2π apart, for which the radial displacement between front and rear contact lines is locally maximal. The capacity of a tube in realisation ω to trap a collar of a given volume is given by the global maximum a + c (ω) among this set, where a ± c (ω) = max 0≤z 0 ≤L a(z 0 ∓ π; ω) − a(z 0 ± π; ω) (III.14) (treating a as an L-periodic function). Likewise if the direction of shear is reversed, the capacity of the tube to trap a collar is determined by a − c (ω).

C. Critical conditions
Equation (III.12) indicates how the shear necessary to support a collar scales with collar volume, supporting the observation (e.g. Figure 3a) that large collars are more resistant to displacement under shear. We therefore assume that the collar contains almost all the fluid in the domain (i.e. V = V 0 ), in order to estimate the largest shear at which a stationary collar can exist in a tube in a given tube realisation ω: This prediction can be tested against the simulation data in Figure 3, for which η = 1: the largest shear was found numerically to be τ c ≈ 0.5045. The small-η asymptotic prediction (III.15) gives the critical shear as τ c ≈ 0.515.
We test these predictions further in Figure 4, using a long section of tube with an irregular shape. For many different values of τ , we solved (III.1) numerically and identified locations at which isolated collars of volume close to V 0 could form stationary equilibria. Different locations were achieved by choosing appropriate initial conditions. Results of simulations show good agreement with the prediction (III.12) (with V = V 0 ), despite the fact that η = 0.5 might be expected to sit outside the range of validity of the small-roughness prediction. The continuous curve in Figure 4(b) arises from evaluating the difference in radial wall location at points an axial distance 2π apart. The segments of this 'finite-difference' curve having positive slope (for which a − z > a + z ) represent stable collar locations both for τ > 0 and when the flow is reversed (τ < 0). The solutions of the algebraic system (III.12) therefore have multiple branches in general; turning points can be identified as saddle-node bifurcation points. In particular, the largest shear at which a collar is recorded coincides with the predicted global maximum given by (III.14, III.15). It is evident from the figure that under quasi-steadily reversing shear, a collar may move smoothly back and forth along a single solution branch if the shear is of small amplitude, jump between branches at larger amplitudes or be swept completely out of the tube at sufficiently large amplitude.
We address the distribution of τ c (ω) over many realisations ω in Sec. V below, after consideration of non-axisymmetric wall roughness.

IV. WEAK NON-AXISYMMETRIC ROUGHNESS
To investigate the effect of non-axisymmetric roughness a(θ, z) on the spatial distribution of annular collars, we have to solve (II.1) numerically, as in Figure 2. However for small roughness and weak shear a perturbation analysis is revealing. We seek steady collar solutions by expanding (II.1) following (III.6), using We assume τ and η are of comparable magnitude as τ → 0, so that the O(η) effects of azimuthal wall roughness appear at O(τ ) in the expansion. The leading-order collar solutions are identical to the axisymmetric case (see (III.9)), with p 0 and p 1 uniform across the collar.
At the following order, Averaging azimuthally using (II.5), integrating with respect to z and imposing zero axial flux at this order gives Thus p 2 is exactly (III.10) and h 2 is given by (III.11) with a replaced by a. The contactline condition h 2 (θ, z ± 0 ) = 0 implies h 2 (z ± 0 ) = 0. Therefore, the collar location in a nonaxisymmetric rough-walled tube must satisfy η a(z 0 − π) − a(z 0 + π) = 12π 3 τ /V, (IV.4) which generalises equation (III.12) to provide a necessary but not sufficient condition for the existence of a stable collar. Likewise we modify (III.14) and (III.15) to define the criticality condition in terms of the azimuthally-averaged tube radius Thus, just as surface tension creates almost axisymmetric structures in a non-axisymmetric tube (Figure 2), we are able to exploit results for an axisymmetric tube to predict outcomes for tubes with two-dimensional roughness.
The present analysis addresses slender collars of length 2π, described by the linearised Young-Laplace equation (II.1b). Larger collars, described by the nonlinear Young-Laplace equation, become shorter as they grow fatter [7,26]; accordingly their stability will be determined by the spatial fields a(z 0 ∓ ; ω) − a(z 0 ± ; ω) for some ≤ 2π.
Having generalised the stability criterion to non-axisymmetric tubes, we now investigate the properties of the maximum asperity a ± c in (IV.5) and evaluate the collar displacement probability (II.9).

V. CRITICAL SHEAR OVER MULTIPLE REALISATIONS
To determine the mean and variance of the critical shear τ c , when sampled over many tube realisations, we can use Monte-Carlo simulation, repeatedly solving (II.1). We can supplement this time-consuming approach by exploiting the algebraic condition (IV.5), valid for small wall roughness, which provides an upper bound on the shear at which stationary collars can be supported within a tube. In addition to applying Monte-Carlo sampling directly and cheaply to (IV.5), we can also approximate the distribution of a c (ω) directly as follows.
The azimuthal averaged random field a(z, ω) is Gaussian with mean zero and covariance (also azimuthally averaged) K 1 increases monotonically with respect to r 1 , from K(r 1 ) ≈ r 1 / √ 2π as r 1 → 0 (an approximation obtained via steepest descents) towards unity for r 1 1. The differenced random field a(z) − a(z + ) (with = 2π for a slender collar) is also Gaussian, with mean zero and Thus the differenced field has variance σ 2 d ≡ K(r 1 )G(0). We normalise the field by defining b(z) = (a(z) − a(z + ))/σ d , which has covariance G(z)/G(0) depending only on r 2 and L.
We then define m = −G zz (0)/G(0), so that 1/ √ m estimates the distance over which the differenced field is correlated. For r 2 L, for example, In this case 1/ √ m ≈ r 2 for r 2 and r 2 / √ 3 for r . Since G is even in in (V.3), it does not distinguish between a + c and a − c , so we shall describe each as a c , which we call the maximum trapping asperity.
We characterise the properties of a c (ω) = σ d max 0≤z<L b(z) in terms of its mean and variance over multiple realisations These are functions of the correlation lengths r 1 and r 2 and the domain length L. It is convenient to exploit the approximation due to Hill et al. [21] of the probability density function P R of max 0≤z<L b as an exponentiated Rayleigh distribution [31], of the form Here N can be interpreted as the effective number of independent local maxima of b in a domain of length L, with each maximum having a Gaussian distribution and with one exceeding all the others in magnitude. For N 1 (i.e. r 2 L), for which terms raised to high powers typically approach either zero or unity, the cumulative distribution function for P R can be approximated to leading order by H(x − (2 log N ) 1/2 ), where H is the Heaviside function (the step sits where P R (x) = 0), implying so that using (II.4)  As r 1 increases, both the mean and variance increase to a plateau (Figure 5a, b), showing how a shorter lengthscale of disorder in the azimuthal direction reduces the capacity of the tube to support stationary collars. The r 1 -dependence arises through the function K 1 , given in (V.2).
As the axial correlation length r 2 increases, the mean critical shear decreases monotonically (figure 5c) whereas its variance increases and then decreases (figure 5d). Clearly an increase in r 2 will reduce axial gradients of the wall, moving the tube towards the smooth case in which a c tends to zero. In contrast, very small r 2 increases the number of potential pinning sites; here the mean depends on log(L/2π) + log(1/r 2 ), indicating a weak dependence on L as r 2 falls. Figure 5 shows good agreement between the mean estimated from the Monte-Carlo simulations (using 10000 samples per data point) from (IV.5) and the approximate probability density function (V.6), and a reasonable prediction of the variance. The simplified estimate (V.8) (assuming N 1 in (V.6)) captures the qualitative form of the mean of a c .
Using (V.6), we can approximate the probability (II.9) that the critical shear is less than a given value τ (corresponding to the probability, over multiple realisations, of the shear τ being sufficient to displace all quasi-steady collars) as Here we assume isotropic covariance (II.4) with correlation length r L; = 2π gives the collar length. Figure 6 shows one example of dependence of the collar displacement probability on the given shear τ . Because determining the critical shear for non-axisymmetric wall perturbation is extremely time consuming, we compare (V.9) with that from Monte-Carlo simulations of the axisymmetric evolution equation equation (III.1), finding respectable agreement between them. The smoothly increasing probability demonstrates the increasing likelihood of eliminating all stationary collars from a tube of a given length as the imposed shear increases; however only in extreme instances can complete removal be guaranteed.
In a very long tube, for example, for which N 1, we expect threshold behaviour with Thus for an isotropic roughness with correlation length r, the threshold shear above which collar removal is almost guaranteed takes the form (1 r L). (V.11) The threshold shear diminishes as r becomes very small (because of the increased disorder in the azimuthal direction) and as r becomes very large (because the tube becomes smoother in the axial direction), with an intermediate maximum where the collar length is comparable in magnitude to the correlation length.

VI. LIQUID PLUGS
Finally, and briefly, a similar argument can be developed for liquid plugs. An isolated plug of radius a 3 0 V in a tube of radius a 0 (1 + ηa(θ, z)) has length a 0 where = V /π + 4 3 to leading order in , where > 2 (so that the plug has hemispherical ends). In the presence of a pressure difference (p − − p + )σ/a 0 , a plug in an axially non-uniform axisymmetric tube can adopt an equilibrium configuration if its hemispherical menisci meet the tube wall with zero contact angle, with the upstream radius being slightly shorter than that downstream. When the wall is non-axisymmetric but weakly rough, we can derive an expression analogous to (IV.5), namely (see Appendix C) We can then exploit directly the estimates of Sec. V, such as the plug displacement probability with N and σ d as given in (V.9). Thus, once again, in a very long tube the displacement probability approaches the step-like form As before, factors inhibiting plug removal include an increase in roughness amplitude η, an increase in tube length and a plug length that is of comparable magnitude to the wall correlation length r.

VII. DISCUSSION
We have investigated the formation and stability of annular collars and plugs in nonaxisymmetric rough-walled tubes, focusing our attention primarily on slender collars, for which a comprehensive analysis is possible. Constrictions in the tube can trap quasi-steady collars, which can then resist displacement under weak shear; likewise, occlusive plugs resist displacement under an imposed pressure gradient. Trapping of collars is reminiscent of nonwetting drops pinned by contact angle hysteresis [5,34], driven drops pinned by surface heterogeneities [37] (and released via sniper bifurcations [41]), wetting drops spreading on rough surfaces [6,38] and capture of fluid from translating droplets by localized asperities [33]. In the present case, surface tension acting through the film's azimuthal curvature promotes localisation of fluid as a collar which fully wets the wall at its effective contact lines. A perturbation analysis for weak wall roughness reveals a simple algebraic condition relating collar location to the imposed shear (IV.4); this appears to be accurate for moderate roughness and is valid for weakly rough non-axisymmetric tubes. From this we deduced a threshold condition for collar trapping in a particular tube realisation (IV.5).
Our primary result is an explicit approximation for the collar displacement probability (V.9) over multiple realisations of a rough-walled tube of finite length, which captures predictions from expensive Monte-Carlo simulations ( Figure 6); it has a straightforward extension for liquid plugs (VI.2) and can be adapted for other scenarios. This approximation was derived in two steps: first, using asymptotic methods to reduce a model framed as a nonlinear partial differential equation to an algebraic system involving some simple geometric features of the irregular wall shape; second, exploitation of a semi-empirical approximation of the distribution of maxima of a Gaussian random process [21]. This approach to uncertainty quantification effectively uses a surrogate model, derived using physical principles, to capture the variability of outcomes (rather than rely on purely numerical approaches), a strategy we have successfully adopted in a related problem [48]. Such low-order descriptions of collar and plug properties are essential in building realistic models of transport in lung airway networks [8,11,36,39], for which stochastic heterogeneity in airway properties is of intrinsic importance.
We focused here primarily on liquid collars, as these are precursors of occlusive liquid plugs in core-annular flows [3,4,26], although our approach is readily adapted to describe plugs directly (Sec. VI). The roughness effects considered here may influence predictions of recent models of liquid plug dynamics in lung airways (e.g. [12,32,42]) under conditions when plug motion is sufficiently slow for the roughness amplitude to be comparable to the trailing film thickness, although numerous additional factors (surfactant, gravity, etc.) will also need to be accounted for. For example our work does not address the case in which the roughness has asperities with axial lengthscales shorter than the film thickness, for which lubrication theory may not be applicable. However our study emphasises that plug or collar dynamics in a real airway (for which geometric properties are known only imperfectly) should sometimes be considered as a stochastic process and clearance of a collar or plug from a real airway can then at best be predicted probabilistically.

Appendix B ASYMPTOTIC STRUCTURE OF A SHEARED COLLAR
We develop an asymptotic approximation of a steady collar under shear by starting in the transition regions, governed by (III.5). We set z = z ± + τ ξ, h = τ 2 H(ξ), q = τ 5 Q, p = P (ξ), so that to leading order Q = 1 2 H 2 + 1 3 H 3 H ξξξ . We assume H → H ± and Q = 1 2 H 2 ± for ξ → ±∞, to match with the adjacent films outside the collar. It is convenient to rescale on these film thicknesses, using H = H ±H (ξ), ξ = H We now consider the solutions of (B.1) where they match the collar, withH 1, satisfying The constants β ± are independent of C ± and reflect translation invariance; C + is determined as a nonlinear eigenvalue (taking the value C + ≈ 2.386) while C − parametrizes the solutions at the rear contact line. We then express (B.2) in terms of the original variables, writing ξ = ξ ± 0 + Z/τ and noting that log |ξ| = log after discarding all terms that are proportional to τ 2 . This motivates the expansion (III. 6) in the outer region.
The first two terms in the outer problem (see (III.9)), can be expanded near each contact line, using Z = z − z 0 ± π, and matched to (B.4) to give Thus the film thickness swept out of the downstream contact line is determined by the collar volume. In the example of Figure 3(b), this prediction gives h shear = τ 2 H + ≈ 0.0272, which is respectably close to the numerically determined value 0.0226 (bearing in mind that the 'small' parameters in this example are τ = 0.5 and η = 1). Likewise, in a strictly steady state, for which H − = H + , we may also set C − = C + to select the structure of the rear contact line.
The τ Z log |Z| and τ Z terms in (B.4) can be matched on to the inner expansion of F (u) in (III.11b) as u → ±π. Crucially, as there is no O(1) contribution to the τ terms in (B.4), the condition h 2 = 0 can then be applied at each contact line to yield (III.12).

Appendix C EQUILIBRIUM LIQUID PLUG IN A ROUGH TUBE
Consider a cylindrical tube with rough wall r = a(θ, z) ≡ 1+ηa(θ, z) in terms of cylindrical coordinates {r, θ, z}; scale lengths on the mean tube radius a 0 and pressure on σ/a 0 . A liquid plug with volume V 0 occupies a domain V within the tube, and is in equilibrium with gas occupying the remainder of the tube. We impose gas pressures p ± downstream (+) and upstream (−) of the plug; the plug's internal pressure p 0 is to be determined. The liquid wets the wall with zero contact angle. The plug's free surfaces S ± have the form r = u ± (θ, z) and touch the wall at z = z ± (θ). Their shape is determined by the Young-Laplace law, contact-line conditions and a volume constraint through ∇ · n ± = −δp ± , ν · n ± = 1 and u ± = a on z = z ± (θ), Here δp ± = p ± − p 0 are pressure differences between the inside and the outside of the free surface, n ± and ν are the unit normals pointing into the gas phase from the free surfaces and the rigid wall respectively, satisfying n ± = (−1, u ±θ /u ± , u ±z )/∆ ± , ν = (−1, a θ /a, a z )/A, (C.2) where ∆ ± (θ, z) = 1 + (u 2 ±θ /u 2 ± ) + u 2 ±z and A(θ, z) = 1 + (a 2 θ /a 2 ) + a 2 z . The mean curvature of the free surfaces is given by When the wall is uniform, a = 1, the free surfaces u ± (z) are hemispherical, satisfying u ± = 1 − (z − z ± ) 2 (z − ≤ z ≤ z − + 1, z + − 1 ≤ z ≤ z + ), δp ± = 2. (C.4) Thus p − = p + and the volume constraint gives the plug length as z + − z − = (V 0 /π) + 4 3 , (C.5) Because z + − z − ≥ 2, we have V 0 ≥ 2π/3. However the plug location remains undetermined.