2.1. 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 $\boldsymbol {u}^*(r^*,z^*,t^*) = u^*\boldsymbol {\hat {r}}+w^*\boldsymbol {\hat {z}}$ and pressure $p^*(r^*,z^*,t^*)$ relative to the gas pressure, which is assumed to be spatially uniform.

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 $\varGamma$. (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\leq r\leq \varPsi _-$) and adjacent to the wall ($\varPsi _+\leq r\leq 1$), with a plug-like region (grey) in between ($\varPsi _-< r<\varPsi _+$).

The liquid is incompressible and we assume inertia is negligible, so conservation of mass and momentum imply

(2.1a,b)\begin{equation} \boldsymbol{\nabla}^*\boldsymbol{\cdot}\boldsymbol{u}^*=\boldsymbol{0},\quad\boldsymbol{\nabla}^*\boldsymbol{\cdot}\boldsymbol{\tau}^* = \boldsymbol{\nabla}^* p^* \quad\text{for}\ R^*\leq r^* \leq a, \end{equation}

where $\boldsymbol {\tau }^*$ is the deviatoric stress tensor. The boundary conditions at the interface are the kinematic condition,

(2.2)\begin{equation} \partial_t^*R^* + w^*\partial^*_zR^* = u^* \quad\mathrm{at}\ r^*=R^*, \end{equation}

and the stress condition,

(2.3)\begin{equation} -p^*\boldsymbol{n} + \boldsymbol{\tau}^*\boldsymbol{\cdot}\boldsymbol{n} = \sigma^*\kappa^*\boldsymbol{n} + \boldsymbol{\nabla}_s^*\sigma^*, \end{equation}

where $\sigma ^*$ is the surface tension, $\kappa ^*$ is the curvature of the interface, $\boldsymbol {n}$ is the unit normal to the interface directed away from the liquid layer and $\boldsymbol {\nabla }_s^*$ is the surface gradient operator. At the tube wall, there is no slip and no penetration,

(2.4)\begin{equation} u^* = w^* = 0 \quad\mathrm{at}\ r^* = a. \end{equation}

The liquid is assumed to be a Bingham fluid with constitutive relation,

(2.5)\begin{equation} \left. \begin{gathered} \tau^*_{ij} = \left(\eta+\frac{\tau_y}{\dot\gamma^*}\right){\dot\gamma}^*_{ij} \quad \mathrm{if}\ \tau^* \geq \tau_y,\\ {\dot\gamma}^*_{ij} = 0 \quad \mathrm{if}\ \tau^*<\tau_y, \end{gathered} \right\} \end{equation}

where $\eta$ is a viscosity, $\tau _y$ is the yield stress, $\boldsymbol {\dot \gamma }^* = \boldsymbol {\nabla }^*\boldsymbol {u}^* + \boldsymbol {\nabla }^*{\boldsymbol {u}^*}^\textrm {T}$ is the shear-rate tensor, and $\tau ^*$ and $\dot \gamma ^*$ are the second invariants of $\boldsymbol {\tau }^*$ and $\boldsymbol {\dot \gamma }^*$, respectively. The second invariant of a tensor $\boldsymbol{\mathsf{T}}$ is defined as $\boldsymbol{\mathsf{T}} = \sqrt {\frac {1}{2}\boldsymbol{\mathsf{T}}_{ij}\boldsymbol{\mathsf{T}}_{ij}}$.

We take the surface tension of the interface to be a linear function of the surfactant concentration,

(2.6)\begin{equation} \sigma^* = \sigma_0 + K(\varGamma_0-\varGamma^*), \end{equation}

where $\sigma _0$ and $\varGamma _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 Reference Stone1990)

(2.7)\begin{equation} \partial^*_t\varGamma^* + \boldsymbol{\nabla}_s^*\boldsymbol{\cdot}(\varGamma^*\boldsymbol{u}^*_s) + \varGamma^*\kappa^*(\boldsymbol{u}^*\boldsymbol{\cdot}\boldsymbol{n})= 0, \end{equation}

where $\boldsymbol {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 Reference Craster and Matar2000; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). At the lateral boundaries of the domain, we impose symmetry boundary conditions,

(2.8)\begin{equation} \partial^*_zR^* = \tau^*_{rz} = w^*=(\boldsymbol{u}^*_s\boldsymbol{\cdot}\boldsymbol{\hat{z}})\varGamma = 0 \quad\mathrm{at}\ z^*=\{0,L^*\}, \end{equation}

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.

2.2. Non-dimensionalisation

To non-dimensionalise (2.1)–(2.7), we introduce the dimensionless quantities,

(2.9)\begin{align} \left. \begin{gathered} (r,z) = \left(\frac{r^*}{a},\frac{z^*}{a}\right), \quad{t} = \frac{\sigma_0}{a\eta}t^*, \quad \boldsymbol{\dot\gamma} = \frac{a\eta}{\sigma_0}\boldsymbol{\dot\gamma}^*, \quad ({u},{w}) = \frac{\eta}{\sigma_0}(u^*,w^*),\quad \varGamma = \frac{\varGamma^*}{\varGamma_0},\\ \boldsymbol{\tau} = \frac{a}{\sigma_0}\boldsymbol{\tau}^*,\quad R = \frac{R^*}{a}, \quad \sigma = \frac{\sigma^*}{\sigma_0}, \quad p = \frac{a}{\sigma_0}p^*, \quad\kappa = a\kappa^*,\quad L = \frac{L^*}{a}. \end{gathered} \right\} \end{align}

The mass and momentum conservation equations (2.1) become

(2.10)$$\begin{gather} 0 = \partial_z{w} + \frac{1}{r}\partial_r(r{u}), \end{gather}$$

(2.11)$$\begin{gather}\partial_r{p} = \frac{1}{r}\partial_r(r\tau_{rr})+\partial_z\tau_{rz}-\frac{\tau_{\theta\theta}}{r}, \end{gather}$$

(2.12)$$\begin{gather}\partial_z{p} = \frac{1}{r}\partial_r(r\tau_{rz}) + \partial_z\tau_{zz}. \end{gather}$$

The wall boundary conditions (2.4) are

(2.13)\begin{equation} {u} = {w} = 0 \quad \text{on} \ r = 1, \end{equation}

the kinematic boundary condition (2.2) is

(2.14)\begin{equation} \partial_{t}R+{w}\partial_zR = {u},\quad\mbox{on}\ r=R, \end{equation}

and the stress boundary condition (2.3) is

(2.15)\begin{equation} -{p}n_i + \tau_{ij}n_j = \sigma\kappa n_i + (\delta_{ij}-n_in_j)\partial_j\sigma, \quad\mbox{on}\ r=R, \end{equation}

where the curvature of the interface is

(2.16)\begin{equation} \kappa = \frac{1}{\sqrt{1+(\partial_zR)^2}}\left[\frac{1}{R}-\frac{\partial_{zz}R}{1+(\partial_zR)^2}\right]. \end{equation}

The constitutive relation (2.5) is

(2.17)\begin{equation} \left. \begin{gathered} \tau_{ij} = \left(1+\frac{\operatorname{\mathcal{B}}}{\dot\gamma}\right){\dot\gamma}_{ij}, \quad \mathrm{if} \ \tau \geq \operatorname{\mathcal{B}},\\ {\dot\gamma}_{ij} = 0 \quad \mathrm{if}\ \tau<\operatorname{\mathcal{B}}, \end{gathered} \right\} \end{equation}

where

(2.18)\begin{equation} {\operatorname{\mathcal{B}}}=\frac{a\tau_y}{\sigma_0} \end{equation}

is a capillary Bingham (or plastocapillarity) number (Jalaal et al. Reference Jalaal, Stoeber and Balmforth2021; van der Kolk et al. Reference van der Kolk, Tieman and Jalaal2023). The equation of state (2.6) becomes

(2.19)\begin{equation} \sigma(\varGamma) = 1 + {\operatorname{\mathcal{M}}}(1-\varGamma), \end{equation}

where

(2.20)\begin{equation} \operatorname{\mathcal{M}} = \frac{K\varGamma_0}{\sigma_0} \end{equation}

is a Marangoni number. The transport equation (2.7) can be expressed in conservative form (Halpern & Frenkel Reference Halpern and Frenkel2003) as

(2.21)\begin{equation} \partial_{t}[R\varGamma\sqrt{1+(\partial_zR)^2}] + \partial_z[{w}_sR\varGamma\sqrt{ 1+(\partial_zR)^2}] = 0, \end{equation}

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, $\mathcal {L}=a/\delta$, where $\delta \ll 1$. We define stretched variables,

(2.22a–f)\begin{equation} \bar{z}=\delta z,\quad \bar{u}=\frac{u}{\delta^2},\quad \bar{w}=\frac{w}{\delta},\quad \bar{t}=\delta^2 {t},\quad \bar{\boldsymbol{\tau}} = \frac{\boldsymbol{\tau}}{\delta},\quad \bar{L} = \delta L, \end{equation}

where the choice of scalings for $z$ and $L$ arise from the long-wave approximation, the scaling for $\boldsymbol {\tau }$ 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. (Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022) 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 $\delta$, 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 Reference Gauglitz and Radke1988; Halpern & Grotberg Reference Halpern and Grotberg1992, Reference Halpern and Grotberg1993; Halpern et al. Reference Halpern, Fujioka and Grotberg2010; Ogrosky Reference Ogrosky2021; Shemilt et al. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022). 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. Reference Camassa, Forest, Lee, Ogrosky and Olander2012; Camassa & Ogrosky Reference Camassa and Ogrosky2015; Ogrosky Reference Ogrosky2021; Shemilt et al. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022), 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 $\delta \ll 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, $|\tau _{rz}|>\operatorname {\mathcal {B}}$ for some $R(z,t)< r<1$, the fluid is fully yielded and the flow is shear-dominated. Where $|\tau _{rz}|\leq \operatorname {\mathcal {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 $\delta$ 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 Reference Walton and Bittleston1991). This general structure is the same as in viscoplastic thin-film flows, as detailed by (Balmforth & Craster Reference Balmforth and Craster1999). 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 $\operatorname {\mathcal {M}}\varGamma _z$ and $p_z$, and their sizes relative to $\operatorname {\mathcal {B}}$. In Appendix A, we show that, in the long-wave limit, the capillary pressure is given by

(2.23)\begin{equation} p = - \kappa\left[1+\operatorname{\mathcal{M}}(1-\varGamma)\right], \end{equation}

where the curvature of the interface, $\kappa$, is defined in (2.16), and that the shear stress is given by

(2.24)\begin{equation} \tau_{rz} = \frac{p_z}{2}\left(r - \frac{R^2}{r}\right) + \frac{R}{r}\mathcal{M}\varGamma_z. \end{equation}

By noting how $\tau _{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\leq r\leq 1$. We define two internal surfaces, $r=\varPsi _-(z,t)$ and $r=\varPsi _+(z,t)$, which separate fully yielded and plug-like regions. The fully yielded regions (where $|\tau _{rz}|>\mathcal {B}$) occupy $R\leq r\leq \varPsi _-$ and $\varPsi _+\leq r\leq 1$, and the plug-like region (where $|\tau _{rz}|\leq \mathcal {B}$) occupies $\varPsi _-\leq r\leq \varPsi _+$. 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).

1. (I) Internal pseudo-plug ($R<\varPsi _-<\varPsi _+<1$). Here, fully yielded regions adjacent to the wall and the interface are separated by a pseudo-plug.

2. (II) Surface pseudo-plug ($\varPsi _-= R<\varPsi _+<1$). Here, the pseudo-plug extends to the interface and the only fully yielded region is adjacent to the wall.

3. (III) Fully yielded ($\varPsi _-=\varPsi _+= R$ or $1=\varPsi _-=\varPsi _+$). In this case, there is no plug-like region and the whole layer is fully yielded.

4. (IV) Near-wall plug ($R<\varPsi _-<1=\varPsi _+$). Here, there is a rigid plug adjacent to the wall and the only fully yielded region is adjacent to the interface.

5. (V) Fully rigid ($\varPsi _-= R$ and $\varPsi _+=1$). Neither yielded region exists and the whole layer is a rigid plug.

Figure 2. (a) 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 $\operatorname {\mathcal {M}}\varGamma _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, $\tilde {w}_s=0$.

Hewitt & Balmforth (Reference Hewitt and Balmforth2012) 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 expressions for $\varPsi _\pm$, from which the type of yielding can be deduced. We find

(2.25)\begin{equation} \varPsi_\pm = \max[R,\min\left(1,\psi_\pm\right)], \end{equation}

where, if $p_z\neq 0$, we have

(2.26) \begin{equation} \psi_{{\pm}} = \left\{\begin{array}{@{}ll} \pm\dfrac{\operatorname{\mathcal{B}}}{|p_z|}+\sqrt{\left(\dfrac{\operatorname{\mathcal{B}}}{p_z}\right)^2+R^2-\dfrac{2R\operatorname{\mathcal{M}}\varGamma_z}{p_z}} & \mbox{if}\ \dfrac{2\operatorname{\mathcal{M}}\varGamma_z}{Rp_z}<1,\\ \dfrac{\operatorname{\mathcal{B}}}{|p_z|}\pm\sqrt{\left(\dfrac{\operatorname{\mathcal{B}}}{p_z}\right)^2+R^2-\dfrac{2R\operatorname{\mathcal{M}}\varGamma_z}{p_z}} & \mbox{if}\ 1+\dfrac{\operatorname{\mathcal{B}}^2}{R^2p_z^2}\geq\dfrac{2\operatorname{\mathcal{M}}\varGamma_z}{Rp_z}\geq1,\\ R & \mbox{if}\ \dfrac{2\operatorname{\mathcal{M}}\varGamma_z}{Rp_z}>1+\dfrac{\operatorname{\mathcal{B}}^2}{R^2p_z^2}, \end{array} \right. \end{equation}

and if $p_z=0$, then $\psi _-=R\operatorname {\mathcal {M}}|\varGamma _z|/\operatorname {\mathcal {B}}$ and $\psi _+=1$. Note that these definitions mean that $\varPsi _\pm$ are continuous everywhere, including at $p_z=0$.

Given values of $R$, $\operatorname {\mathcal {B}}$, $p_z$ and $\operatorname {\mathcal {M}}\varGamma _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/\operatorname {\mathcal {B}}$, and on the size of the Marangoni stress relative to the yield stress, via $\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {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/\operatorname {\mathcal {B}},\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {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 $\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {B}}$ is then decreased, so that we cross the line $\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {B}}=1$. Then the shear stress at the interface has dropped below $\operatorname {\mathcal {B}}$ so there can no longer be a fully yielded region at the interface, and the layer then exhibits yielding of type II. Decreasing $\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {B}}$ further so that $\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {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/\operatorname {\mathcal {B}}$, we remain in yielding type III until we reach the line $(1-R^2)p_z+2R\operatorname {\mathcal {M}}\varGamma _z=-2\operatorname {\mathcal {B}}$. 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/\operatorname {\mathcal {B}}$ yet further, so that we cross the line $(1-R^2)p_z+2R\operatorname {\mathcal {M}}\varGamma _z=2\operatorname {\mathcal {B}}$, now $p_z/\operatorname {\mathcal {B}}$ 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/\operatorname {\mathcal {B}}$ and $\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {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 $(p_z/\operatorname {\mathcal {B}},\operatorname {\mathcal {M}}\varGamma _z/\operatorname {\mathcal {B}})=(\pm 2/(1-R),\pm 1)$. 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 $\varPsi _\pm$ 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

(2.27)\begin{equation} R_t=\frac{1}{R}Q_z, \end{equation}

where the axial volume flux is given by

(2.28) \begin{equation} Q = \left\{\begin{array}{@{}ll} -\dfrac{p_z}{16}F_1-\dfrac{1}{4}R\operatorname{\mathcal{M}}\varGamma_zF_2-\dfrac{\operatorname{\mathcal{B}}}{6}\operatorname{sgn}(p_z)(F_3+F_4) & \mbox{if}\ \dfrac{2\operatorname{\mathcal{M}}\varGamma_z}{Rp_z}<1,\\ -\dfrac{p_z}{16}F_1-\dfrac{1}{4}R\operatorname{\mathcal{M}}\varGamma_zF_2-\dfrac{\operatorname{\mathcal{B}}}{6}\operatorname{sgn}(p_z)(F_3-F_4) & \mbox{if}\ \dfrac{2\operatorname{\mathcal{M}}\varGamma_z}{Rp_z}\geq 1,\\ -\dfrac{1}{4}R\operatorname{\mathcal{M}}\varGamma_zF_2 + \dfrac{\operatorname{\mathcal{B}}}{6}\operatorname{sgn}(\varGamma_z)F_4 & \mbox{if}\ p_z=0, \end{array} \right. \end{equation}

with

(2.29a)$$\begin{gather} F_1 = \varPsi_-^4-4R^2\varPsi_-^2+R^4\left[3-4\log\left(\frac{R\varPsi_+}{\varPsi_-}\right)\right] -\varPsi_+^4 + 4R^2\varPsi_+^2-4R^2+1, \end{gather}$$

(2.29b)$$\begin{gather}F_2 = \varPsi_-^2-R^2-\varPsi_+^2+1+2R^2\log\left(\frac{R\varPsi_+}{\varPsi_-}\right), \end{gather}$$

(2.29c)$$\begin{gather}F_3 = (\varPsi_+ - 1)( - 3R^2+1+\varPsi_+ + \varPsi_+^2),\quad F_4=\varPsi_-^3-3R^2\varPsi_- + 2R^3. \end{gather}$$

If $\operatorname {\mathcal {M}}=0$, or equivalently if there is no surfactant, $\varGamma =0$, then $\varPsi _-=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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022). The surfactant transport equation is

(2.30)\begin{equation} (R\varGamma)_t + (w_sR\varGamma)_z=0,\end{equation}

where the surface velocity is

(2.31) \begin{equation} w_s = \left\{\begin{array}{@{}ll} \dfrac{1}{4}p_zG_1+R\operatorname{\mathcal{M}}\varGamma_zG_2+\operatorname{\mathcal{B}}\operatorname{sgn}(p_z)(G_3+G_4) & \mbox{if}\ \dfrac{2\operatorname{\mathcal{M}}\varGamma_z}{Rp_z}<1,\\ \dfrac{1}{4}p_zG_1+R\operatorname{\mathcal{M}}\varGamma_zG_2+\operatorname{\mathcal{B}}\operatorname{sgn}(p_z)(G_3-G_4) & \mbox{if}\\dfrac{2\operatorname{\mathcal{M}}\varGamma_z}{Rp_z}\geq1,\\ R\operatorname{\mathcal{M}}\varGamma_zG_2 - \operatorname{\mathcal{B}}\operatorname{sgn}(\varGamma_z)G_4 & \mbox{if}\ p_z=0, \end{array}\right. \end{equation}

with

(2.32a)\begin{gather} G_1 = R^2+\varPsi_+^2-\varPsi_-^2-1-2R^2\log\left(\frac{R\varPsi_+}{\varPsi_-}\right), \end{gather}

(2.32b–d)\begin{gather} G_2 = \log\left(\frac{R\varPsi_+}{\varPsi_-}\right),\quad G_3 = 1-\varPsi_+,\quad G_4 = R-\varPsi_-. \end{gather}

Equations (2.27)–(2.32) are solved subject to the boundary conditions

(2.33)\begin{equation} R_z = Q = w_s\varGamma = 0, \quad\mbox{at}\ z=\{0,L\}, \end{equation}

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)\begin{equation} \tau_w = \frac{p_z}{2}(1-R^2) + R\operatorname{\mathcal{M}}\varGamma_z. \end{equation}

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|\ll 1$, but no longer require that $\delta \ll 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

(2.35)\begin{equation} R(z,t) = 1 - \epsilon H(z,t), \end{equation}

where $\epsilon \ll 1$ and $H$ is the scaled layer thickness. Similarly, we define the boundaries between fully yielded and plug-like regions,

(2.36)\begin{equation} \varPsi_\pm = 1 - \epsilon Y_\mp. \end{equation}

Since $\varPsi _+\geq \varPsi _-$, the definitions (2.36) mean $Y_+\geq Y_-$. Other relevant variables are rescaled by defining

(2.37)\begin{equation} \tilde{p} = \frac{1+p}{\epsilon},\quad\tilde{\kappa} = \frac{\kappa-1}{\epsilon},\quad \tilde{t}=\epsilon^3t, \quad {\tilde{\boldsymbol{\tau}} = \frac{\boldsymbol{\tau}}{\epsilon^2}}, \end{equation}

and the scaled capillary Bingham and Marangoni numbers are given respectively by

(2.38)\begin{equation} {B} = \frac{\operatorname{\mathcal{B}}}{\epsilon^2} = \frac{a\tau_y}{\sigma_0\epsilon^2}\quad\mathrm{and}\quad M = \frac{\operatorname{\mathcal{M}}}{\epsilon^2} = \frac{K\varGamma_0}{\sigma_0\epsilon^2}. \end{equation}

We then insert (2.35)–(2.38) into the long-wave equations (2.23)–(2.33) and truncate at leading order in $\epsilon$. From (2.23), the thin-film capillary pressure gradient is given by

(2.39)\begin{equation} \tilde{p}_z = - \tilde\kappa_z = - H_z-H_{zzz} \end{equation}

to leading order in $\epsilon$. (Unlike in the thick-film case (2.23), surfactant has no effect on the mean surface tension at leading order in the thin-film limit.) Note that with the scalings (2.37) and (2.38), $2\operatorname {\mathcal {M}}\varGamma _z/R{p}_z=O(\epsilon )$, so in the thin-film theory, we always have $2\operatorname {\mathcal {M}}\varGamma _z/R{p}_z<1$. This simplifies (2.26) somewhat, giving $Y_\pm =\max [0,\min (H,\mathcal {Y}_\pm )]$ where

(2.40)\begin{equation} \mathcal{Y}_\pm = H + \frac{{M}\varGamma_z}{p_z}\pm\frac{B}{|p_z|}, \end{equation}

assuming for now that $p_z\neq 0$. As in § 2.3, we can present criteria for each of the five yielding types to occur based on the values of $Y_\pm$ (as in figure 2a): I, internal pseudo-plug ($0< Y_-< Y_+< H$); II, surface pseudo-plug ($0< Y_-< H=Y_+$); III, fully yielded ($Y_\pm =0$ or $Y_\pm =H$); IV, near-wall plug ($Y_-=0< Y_+< H$); V, fully rigid ($Y_-=0$ and $Y_+=H$). Figure 2(c) illustrates where in $(H\tilde {p}_z/B,M\varGamma _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.

The evolution equation (2.27), to leading order in $\epsilon$, becomes

(2.41)\begin{equation} H_{\tilde{t}} + q_z = 0, \end{equation}

where the scaled axial volume flux is

(2.42)$$\begin{gather} q = - \frac{1}{3}\tilde{p}_z[H^3+(H-Y_+)^3-(H-Y_-)^3]-\frac{1}{2}M\varGamma_z[H^2-(H-Y_-)^2+(H-Y_+)^2]\nonumber\\ +\frac{1}{2}B\operatorname{sgn}(\tilde{p}_z)[H^2-(H-Y_-)^2-(H-Y_+)^2], \end{gather}$$

if $\tilde {p}_z\neq 0$. When $\tilde {p}_z=0$, the flux is $q=-\tfrac {1}{2}\operatorname {sgn}(\varGamma _z)H^2(|M\varGamma _z|-B)$ if $|M\varGamma _z|>B$, and $q=0$ if $|M\varGamma _z|\leq B$. As can be seen in figure 2(c), when $\tilde {p}_z=0$, the only type of yielding possible is type III, which occurs if $|M\varGamma _z|>B$; the layer is rigid (type V) if $|M\varGamma _z|\leq B$.

The surfactant transport equation (2.30), at leading order, becomes

(2.43)\begin{equation} \varGamma_{\tilde{t}} + [\tilde{w}_s\varGamma]_z = 0, \end{equation}

where the scaled surface velocity is

(2.44)$$\begin{gather} \tilde{w}_s = - \frac{1}{2}\tilde{p}_z[H^2+(H-Y_+)^2-(H-Y_-)^2] - M\varGamma_z(H+Y_- - Y_+)\nonumber\\ -B\operatorname{sgn}(\tilde{p}_z)(H-Y_- - Y_+), \end{gather}$$

if $\tilde {p}_z\neq 0$. When $\tilde {p}_z=0$, the surface velocity is $\tilde {w}_s = -\operatorname {sgn}(\varGamma _z)H(|M\varGamma _z|-B)$ if $|M\varGamma _z|>B$ and $\tilde {w}_s=0$ if $|M\varGamma _z|\leq B$. From (2.40) and (2.44), we can deduce that if $H\tilde {p}_z=-2M\varGamma _z$, then $Y_-=H-Y_+$ and $\tilde {w}_s=0$, meaning that the interface of the layer is immobilised. 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)\begin{equation} H_z = q = \varGamma \tilde{w}_s = 0, \quad\mathrm{at} \ z= \{0,L\}. \end{equation}

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)\begin{equation} \tilde{\tau}_w = H\tilde{p}_z + M\varGamma_z. \end{equation}

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

(2.47)\begin{equation} H(z,t=0) = 1 - A \cos\left(\frac{{\rm \pi} z}{L}\right),\quad \varGamma(z,t=0) = 1, \end{equation}

and for the long-wave equations, the equivalent conditions are

(2.48)\begin{equation} R(z,t=0) = \sqrt{(1-\epsilon)^2-\epsilon^2 A^2/2}-\epsilon A\cos\left(\frac{{\rm \pi} z}{L}\right), \quad \varGamma(z,t=0) = 1, \end{equation}

where $A$ is the amplitude of the initial perturbation and $\epsilon$ 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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022), 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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022) and other related studies (Gauglitz & Radke Reference Gauglitz and Radke1988; Halpern et al. Reference Halpern, Fujioka and Grotberg2010), we solve the equations in a domain of length ${L=\sqrt {2}{\rm \pi} }$. This is the domain length that can accommodate the most unstable mode in the Newtonian linear stability analysis (Hammond Reference Hammond1983). 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 (Reference Jalaal2016), which was used for the surfactant-free problem (Shemilt et al. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022) and has also been used for several other viscoplastic thin-film problems (Balmforth et al. Reference Balmforth, Ghadge and Myers2007; Jalaal et al. Reference Jalaal, Stoeber and Balmforth2021). We define $\hat {Y}_\pm =\max (Y_{min},Y_\pm )$ and $\hat {\varPsi }_\pm =\min (1-Y_{min},\varPsi _\pm )$ where $Y_{min}$ is a small parameter, and replace $Y_\pm$ and $\varPsi _\pm$ with $\hat {Y}_\pm$ and $\hat {\varPsi }_\pm$, 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\leq z\leq 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.

3.1. Time evolution of a surfactant-laden thin film

Figure 3(a) shows snapshots from a sample numerical solution of the thin-film system. At $\tilde {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 $\tilde {t}=5$. In this very-early-time period, $\varGamma _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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022). 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.

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 $\tilde {t}=\{0,5,80,100,130,250,500,2000\}$. At each $\tilde {t}$, there are three panels: the top panel shows the layer evolving, with $Y_-$ (cyan) and $Y_+$ (red), and the thin-film axial velocity $\tilde {w}$ represented by the colour map; the middle panel shows plots of $\tilde {p}_z$ (magenta), $\varGamma$ (green) and $10\varGamma _z$ (blue); the bottom panel shows the solution in $(H\tilde {p}_z/B$,$M\varGamma _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 _zH$ from the same simulation (solid) compared with the evolution of $\max _zH$ 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 _zY_-$ (solid red), $\max _z(H-Y_+)$ (solid blue) and $\max _z(|\tilde {\tau }_w|-B)$ (solid black) for the simulation in panel (a). Also shown are plots of $\max _zY_ - $ (dash-dotted red) and $\max _z(|\tilde {\tau }_w|-B)$ (dash-dotted black) for the surfactant-free viscoplastic simulation $(B=0.04,M=0)$.

Between $\tilde {t}=5$ and $\tilde {t}=80$ (figure 3a), significant gradients in surfactant concentration develop, and at $\tilde {t}\approx 80$, a fully yielded region at the interface, where $M\varGamma _z$ exceeds $B$, appears near the centre of the domain. This yielded region grows and has extended across the whole interface by $\tilde {t}=500$. During the intermediate period, $80<\tilde {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 $\varGamma _z$. At $\tilde {t}=130$, the right-travelling wave in $\varGamma _z$ near $z=4$ resembles a discontinuous shock wave. The sudden decrease in $\varGamma _z$ across the shock results in a rise in $Y_-$ just ahead of the shock where the Marangoni force is weaker. Similarly, at $\tilde {t}=250$, a shock-like discontinuity has developed in $\varGamma _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 $\varGamma _z$, we can use (2.43) to derive a Rankine–Hugoniot condition (see Appendix B), which provides the relation

(3.1)\begin{equation} u_s = \frac{[(\tilde{w}_s\varGamma)_z]_-^+}{[\varGamma_z]_-^+}\end{equation}

for the shock propagation speed, $u_s$, in terms of the sizes of jumps in $\varGamma _z$ and $(\tilde {w}_s\varGamma )_z$ across the discontinuity. We have already noted how the jump in $\varGamma _z$ across the shock affects the yielding behaviour on either side, but (3.1) also shows that the yielding behaviour, via $\tilde {w}_s$, 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 $\tilde {t}=80$, $\tilde {t}=100$ and $\tilde {t}=130$, at both side boundaries, the solution in $(H\tilde {p}_z$/B,$M\varGamma _z/B)$-space approaches $(H\tilde {p}_z/B,M\varGamma _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_-\rightarrow 0$ as $z\rightarrow 0$ or $z\rightarrow 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\rightarrow 0$ or $z\rightarrow L$ during this period is given by $\lim _{z\rightarrow 0,L}(H-M\varGamma _{zz}/\tilde {p}_{zz})$, which can be seen to take a positive value in figure 3(a) for $80\leq \tilde {t}\leq 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 layer while getting gradually smaller, indicating that the layer is rigidifying. Both $Y_-$ and $H-Y_+$ tend towards zero at a rate proportional to $1/t$ (figure 3c). It can also be seen that in $(H\tilde {p}_z/B,M\varGamma _z/B)$-space (figure 3a, $\tilde {t}=500,2000$), the whole solution is close to the line $M\varGamma _z+\tfrac {1}{2}H\tilde {p}_z=0$. This indicates that $H-Y_+\approx Y_-$, and so $\tilde {w}_s\approx 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\tilde {p}_z/B,M\varGamma _z/B)=(-2,1)$ (figure 3a, $\tilde {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\rightarrow H_0(z;B)$ as $\tilde {t}\rightarrow \infty$, which satisfies

(3.2a)\begin{equation} H_0(H_{0,z}+H_{0,zzz}) = 2B,\end{equation}

whilst $M\varGamma \rightarrow M\varGamma _0$ as $\tilde {t}\rightarrow \infty$ where $\varGamma _0$ has a linear profile with slope $B/M$ and mass $L$ over $0< z< L$,

(3.2b)\begin{equation} M\varGamma_0 = M - \frac{BL}{2} + Bz. \end{equation}

For the solution (3.2b) to be valid, it must have $\varGamma _0\geq 0$ everywhere, or specifically $2M\geq BL$. Indeed, we observe that in simulations with $2M< BL$, $H$ and $\varGamma$ do not approach $H_0$ and $\varGamma _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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022), 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(|\tilde {\tau }_w|)$, the maximum value of the shear stress exerted on the tube wall, showing that it peaks around $\tilde {t}\approx 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.

3.2. Late-time dynamics of a thin film

To predict the late-time behaviour of the layer, we propose an asymptotic solution for $\tilde {t}\gg 1$. Inspired by observations from numerical simulations such as figure 3, we make expansions of the form

(3.3)\begin{equation} \left. \begin{gathered} H = H_0 + \frac{H_1}{Bt}+\cdots,\quad \varGamma = \varGamma_0 + \frac{\varGamma_1}{Bt}+\cdots,\quad Y_+ = H_0 + \frac{Y_{+,1}}{Bt}+\cdots,\\ Y_- = \frac{Y_{-,1}}{Bt}+\cdots, \quad \partial_zp = \partial_zp_0 + \frac{\partial_zp_1}{Bt}+\cdots = - \frac{2B}{H_0} - \frac{(\partial_z+\partial_z^3)H_1}{Bt} + \cdots, \end{gathered} \right\} \end{equation}

where $\varGamma _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\geq 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. Substituting (3.3) into (2.40)–(2.42) gives

(3.4a–c)\begin{align} H_1 = \partial_z(Y_{-,1}^2),\quad Y_{+,1} = H_1 - \frac{H_0}{2B}M\partial_z\varGamma_{1},\quad Y_{-,1} = Y_{+,1} + \frac{H_0^2}{2B}(\partial_z+\partial_z^3)H_1. \end{align}

Inserting (3.3) into the surface velocity (2.44) gives, to leading order,

(3.5)\begin{equation} \tilde{w}_s \sim \frac{1}{B^2t^2}W_1, \quad \mathrm{where}\ W_1 = \frac{BY_{-,1}^2}{H_0}-\frac{H_0}{4B}(M\partial_z\varGamma_{1})^2, \end{equation}

and the transport equation (2.43) gives

(3.6)\begin{equation} B\varGamma_1 = \partial_z\left(\varGamma_0W_1\right). \end{equation}

The boundary conditions (2.45) imply

(3.7)\begin{equation} \partial_zH_0 = \partial_zH_1 = Y_{-,1} = W_1 = 0, \quad\mathrm{at}\ z= \{0,L\}, \end{equation}

and mass conservation implies

(3.8)\begin{equation} \int_0^LH_0\,\mathrm{d}z = L,\quad \int_0^LH_1\,\mathrm{d}z = 0, \quad\int_0^L\varGamma_1\,\mathrm{d}z=0. \end{equation}

We solve (3.2a) and (3.4)–(3.6) numerically subject to (3.7) and (3.8).

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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022) 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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022). 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. (a) Plot of $\max _zH_0$ for all the solutions to (3.2a) (solid), which are static, marginally yielded solutions to the thin-film equations with surfactant. For comparison, the equivalent plot for the static solutions for the surfactant-free problem, which were computed by Shemilt et al. (Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022), is also shown (dotted). The four coloured markers in panel (a) correspond to the example solutions shown in panel (b). All solutions in panel (b) have $B=0.05$ but the dotted ones are the solutions for the surfactant-free problem. (c) $O(1/t)$ terms in the late-time expansion (3.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 $\varGamma _1$ (solid blue). These are compared with the corresponding quantities from a numerical solution of the thin-film equations (2.40)–(2.45) at $\tilde {t}=10^4$, specifically $[H(z,\tilde {t}=10^4)-H_0(z)]B\tilde {t}$ (dashed black), $Y_-(z,\tilde {t}=10^4)B\tilde {t}$ (dashed magenta), $[Y_+(z,\tilde {t}=10^4)-H(z,\tilde {t}=10^4)]B\tilde {t}$ (dashed red) and $[\varGamma (z,\tilde {t}=10^4)-\varGamma _0(z)]B\tilde {t}$ (dashed cyan) where $\varGamma _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\}$.

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 $\tilde {t}=10^4$. Also, figure 4(c) suggests that $Y_{-,1}\approx 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, $\tilde {w}_s$, 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.3. 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 _zH(z,\tilde {t}=10^4)$, which have several characteristics of note.

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\leq \max _zH\leq 2.8$. The red lines are $2M=BL$ (solid), $B=B_m(A=0.2)\approx 0.0289$ (dashed) and $B=2B_m(A=0.2)\approx 0.0578$ (dash-dotted).

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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022), 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. (Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022) 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)\approx 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}\approx 2B_m\approx 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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022), 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\varGamma _z\approx 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\varGamma _z\approx 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\gg 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,

(3.9)\begin{equation} H_{\tilde{t}} + \left[\frac{1}{6}\tilde{p}_z\tilde{Y}^2(2\tilde{Y}-3H)\right]_z=0, \quad\mbox{where}\ \tilde{Y} = \max\left(0,\frac{1}{2}H-\frac{B}{|\tilde{p}_z|}\right). \end{equation}

If we replace $B$ with $\tfrac {1}{2}B$ and $\tilde {t}$ with $\tfrac {1}{4}\tilde {t}$ in (3.9), then we exactly recover the evolution equation for the surfactant-free problem ((2.27) of Shemilt et al. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022). 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. Reference Otis, Johnson, Pedley and Kamm1993), 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\gg 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\geq 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\geq 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< B L$. This means that the layer cannot enter the late-time regime described in § 3.2, as then (3.2b) would require $\varGamma$ 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=B L$. 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\rightarrow 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 $\tilde {p}_z$ 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(\epsilon ^2)$) reduction in surface tension relative to its mean value. In this thin-film limit, surfactant operates solely through powerful Marangoni effects.

Figure 6. Snapshots from a thin-film simulation with $B=0.04$, $M=0.08$, $A=0.2$, at $\tilde {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, $\tilde {w}$, represented by the colour map; the middle row shows $\varGamma$ (green), $\varGamma _z$ (blue) and $\tilde {p}_z$ (magenta); and the lower row shows shows the solution in $(H\tilde {p}_z/B,M\varGamma _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. (Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022) 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 _zH(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. 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 _zH(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).

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.

4.1. Example evolution showing plug formation

Figure 8(a) shows a numerical solution of the long-wave equations. The layer thickness is $\epsilon =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 $\tilde {t}\approx 268$. Around $267\lesssim \tilde {t}\lesssim 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 Reference Ogrosky2021). 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) Numerical solution of the long-wave equations (2.25)–(2.33) with $A=0.2$, $\operatorname {\mathcal {B}}=0.001$, $\mathcal {M}=0.02$ and $\epsilon =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), $\varGamma _z$ (blue) and $\varGamma$ (green). (b) Time evolution of $\max _zH$ for the same simulation, compared with $\max _zH$ from simulations with $(\operatorname {\mathcal {B}},\operatorname {\mathcal {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.

Figure 8(a) shows that by $\tilde {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. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022). 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.

4.2. 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 $\operatorname {\mathcal {M}}\gg 1$,

(4.1)\begin{align} R = \mathcal{R}_0 + \frac{1}{\operatorname{\mathcal{M}}}\mathcal{R}_1+\cdots,\quad \varGamma = \mathcal{G}_0 + \frac{1}{\operatorname{\mathcal{M}}}\mathcal{G}_1+\cdots,\quad w_s = \mathcal{W}_0 + \frac{1}{\operatorname{\mathcal{M}}}\mathcal{W}_1+\cdots. \end{align}

Since the surface stress must be finite, $|\partial _z\sigma |=|\operatorname {\mathcal {M}}\partial _z\varGamma |<\infty$, we must have $\partial _z\mathcal {G}_0=0$. Conservation of the total amount of surfactant in the long-wave theory then implies

(4.2)\begin{equation} \mathcal{G}_0(t) = \frac{\int_0^L\mathcal{R}_0|_{t=0}\,\mathrm{d}z}{\int_0^L\mathcal{R}_0\,\mathrm{d}z}. \end{equation}

The surfactant transport equation (2.30) at leading order in $\operatorname {\mathcal {M}}$, after rearranging and using (4.2), gives

(4.3)\begin{equation} \mathcal{W}_0(z,t) = \frac{1}{\mathcal{R}_0}\int_0^z\mathcal{I}(\zeta,t)\,\mathrm{d}\zeta,\quad\mbox{where}\ \mathcal{I} = \frac{\mathcal{R}_0\int_0^L\partial_t\mathcal{R}_{0}\,\mathrm{d}z-\partial_t\mathcal{R}_{0}\int_0^L\mathcal{R}_0\,\mathrm{d}z}{\int_0^L\mathcal{R}_0\,\mathrm{d}z}. \end{equation}

Therefore, unlike for the thin-film case (Appendix C), the surface velocity is not zero for large $\operatorname {\mathcal {M}}$, in general. From (4.3), the surface velocity can be deduced from the difference between the local time derivative of the film thickness, $\partial _t\mathcal {R}_0$, and the globally averaged rate of change of the thickness, via the integral of $\partial _t\mathcal {R}_0$. Non-zero $\mathcal {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 $\sigma '(\varGamma )\neq 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 $\sigma '(\varGamma )\neq 0$.

Figure 9 compares a solution of the full long-wave system (2.27)–(2.33) with $\operatorname {\mathcal {M}}=10$ with the results (4.2) and (4.3). It shows good agreement between $\varGamma$ and $w_s$ from the numerical solutions and the leading-order approximations (4.2) and (4.3), despite the value of $\mathcal {M}$ not being extremely large. The spatial gradients in $\varGamma$ (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-$\operatorname {\mathcal {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.

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

The results from figure 9 suggest that the large-$\operatorname {\mathcal {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 $\operatorname {\mathcal {M}}$ in simulations while retaining accuracy. Figure 9(b) suggests that when the simulation is stopped, the increase in $\varGamma$ from its original value was approximately $3\,\%$. We have found that this is approximately independent of $\operatorname {\mathcal {M}}$, so for values of $\operatorname {\mathcal {M}}$ close to or larger than $\operatorname {\mathcal {M}}\approx 30$, we expect the theory to break down since $\sigma$ will approach zero as a plug forms. For this reason, $\operatorname {\mathcal {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.

4.3. Critical layer thickness for plug formation

Figure 10(a) compares the computed critical layer thickness for plug formation to occur, $\epsilon _{crit}$, when surfactant is present with its value when the interface is clean. Both with and without surfactant, increasing $\operatorname {\mathcal {B}}$ induces an increase in $\epsilon _{crit}$. The increase is small when $\operatorname {\mathcal {B}}$ is small, but is significant at larger $\operatorname {\mathcal {B}}$. When surfactant is present, the increase in $\epsilon _{crit}$ is initiated at $\operatorname {\mathcal {B}}\approx 10^{-3}$, compared with at $\operatorname {\mathcal {B}}\approx 2\times 10^{-3}$ when the interface is clean. Beyond these values of $\operatorname {\mathcal {B}}$, the increase is amplified by the presence of surfactant, such that the value of $\epsilon _{crit}$ is approximately $30\,\%$ larger than when the interface is clean.

Figure 10. (a) Critical layer thickness, $\epsilon _{crit}$, as a function of $\mathcal {B}$, for a surfactant-free layer ($\mathcal {M}=0$) and a layer with surfactant ($\mathcal {M}=0.4$). All simulations have $A=0.25$. The value of $\epsilon _{crit}$ computed is such that a simulation with $\epsilon =\epsilon _{crit}+0.001$ forms a plug before $\tilde {t}=10^4$ and a simulation with $\epsilon =\epsilon _{crit}-0.001$ has not formed a plug by $\tilde {t}=10^4$. (b) Data from long-wave simulations at various values of $\epsilon$ and $\operatorname {\mathcal {M}}$, with $A=0.25$ and $\operatorname {\mathcal {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$.

Figure 10(b) illustrates the effect on the dynamics of varying $\operatorname {\mathcal {M}}$ and $\epsilon$. 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 $\epsilon _{crit}$, so the dependence of the critical layer thickness on $\operatorname {\mathcal {M}}$ can be seen. At low $\operatorname {\mathcal {M}}$, the critical thickness for the surfactant-free problem is recovered, and as $\operatorname {\mathcal {M}}$ is increased, $\epsilon _{crit}$ increases to a significantly higher value for large $\operatorname {\mathcal {M}}$. Therefore, for this value of $\operatorname {\mathcal {B}}$, there is a range of layer thicknesses $0.14\lesssim \epsilon \lesssim 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}<\operatorname {\mathcal {M}}<10^{-1}$, over which $\epsilon _{crit}$ changes appreciably in figure 10(b), corresponds to a maximum relative reduction of surface tension of between $1\,\%$ and $10\,\%$. Figure 10(b) also quantifies how the plugging time increases as $\operatorname {\mathcal {M}}$ is increased. We know increasing $\operatorname {\mathcal {B}}$ extends the plugging time (Shemilt et al. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022), 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 Reference Halpern and Grotberg1993). Instead, figure 10(b) highlights how the increase in $\epsilon _{crit}$ due to rigidification and stabilisation of the layer by yield stress is amplified by the presence of sufficiently strong surfactant.

The value of $\epsilon _{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 $\operatorname {\mathcal {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, $\tilde {t}=10^4$, will have some effect on the computed value of $\epsilon _{crit}$, but figure 10(b) shows that the plugging time increases rapidly close to the computed value of $\epsilon _{crit}$, suggesting that there is good convergence to the true value of $\epsilon _{crit}$.