2. Mathematical model

2.1. Long-wave equations

We model the airway as a straight cylindrical tube of radius $R$ , lined by an annular film of mucus, which is assumed to be a Newtonian fluid with viscosity ${\mu} _m$ and density $\rho _m$ . Air of viscosity ${\mu} _a$ and density $\rho _a$ flows through the core. The mucus film interacts with air through an interface that is endowed with interfacial tension $\gamma$ . We assume axisymmetry and adopt a cylindrical coordinate system (see figure 1 b). The axial length scale $\varLambda$ that characterises axial variations is assumed to be much larger than the radius $R$ so that $\epsilon = R/\varLambda$ is a small parameter. We can then treat the two-phase flow in the thin-film limit and obtain simplified long-wave equations. Following Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015), we scale the continuity and Navier–Stokes equations, and retain terms up to $\mathcal{O}(\epsilon ^2)$ to obtain the following non-dimensional, long-wave, boundary-layer equations (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012) for pressure, $p_i$ , and the radial and axial velocities, $v_i$ and $u_i$ :

(2.1) \begin{align}\frac {1}{r}{\partial _r}(r v_i) + {\partial _z u_i} &= 0,\end{align}

(2.2) \begin{align} \epsilon {\textit{Re}}_i\left (\partial _t u_i + v_i\partial _r u_i+ u_i\partial _z u_i\right ) &= -\partial _z \left (p_i |_{d}\right ) -\epsilon ^2 \partial _z \left (\partial _z u_i |_{d}\right ) +\frac {1}{r}\partial _r \left (r\partial _r u_i\right ) \nonumber \\ &\quad +\, 2\epsilon ^2 \partial _{zz} u_i+ G_{i}\sin( \omega t), \end{align}

where all variables are non-dimensional and $\partial _t$ represents the partial derivative with respect to time $t$ , and so on for $\partial _r$ and $\partial _z$ . The interface between the air and mucus (denoted by subscripts $i = a, \, m$ ) is located at $r = d(z,t)$ . The two terms evaluated at the interface, i.e. those involving $p_i |_{d}$ and $u_i |_{d}$ , were obtained by substituting for $\partial _z p_i$ the expression obtained after integrating the $\mathcal{O}(\epsilon ^2)$ radial momentum equation from $r$ to $d$ .

The following characteristic scales (decorated by ^*) have been used for non-dimensionalisation:

(2.3) \begin{equation} r^* = {R},\quad z^* ={\varLambda },\quad t^* = {R/U},\quad, d^* = {R},\, u^*_i ={U}, \quad v^*_i =\epsilon {U}, \quad p_i^* = \frac {{\mu} _i U}{\epsilon R}, \end{equation}

where $U$ , the common velocity scale for both phases, will be chosen shortly.

The last term of (2.2) arises from the imposition of an oscillatory body force (or a background, uniform, axial pressure-gradient) in the air phase ( $G_m = 0$ ), of frequency $f_b =\omega /2 \pi$ and dimensional amplitude $A = G_a {\mu} _a U/R^2$ , to mimic respiratory airflow. We now select the velocity scale $U$ based on the airflow, by setting $G_a = 1$ so that $U = A R^2/{\mu} _a$ . Then, the non-dimensional Reynolds numbers, ${\textit{Re}}_i = \rho _i \textit{UR}/{{\mu} _i}$ , become ${\textit{Re}}_a= \rho _a A R^3/ {\mu} _a^2$ and ${\textit{Re}}_m = {\textit{Re}}_a \varPi _{{\mu} }/\varPi _{\rho }$ , where $\varPi _{{\mu} } = {\mu} _a/{\mu} _m$ and $\varPi _{\rho } = \rho _a/\rho _m$ are the viscosity and density ratios.

Next, we list the conditions at the boundaries, starting with the interface, $r = d(z,t)$ , where we require the velocities to be continuous,

(2.4) \begin{align} u_a = u_m , \quad v_a = v_m, \quad \textit{at}\quad r = d, \end{align}

and apply the $\mathcal{O}(\epsilon ^2)$ balances of tangential and normal stresses:

(2.5) \begin{align}\partial _r u_m-\varPi _{\mu} \partial _r u_a &= \left [2\epsilon ^2 \partial _z d\left (\partial _z u_m- \partial _r v_m\right )-\epsilon ^2 \partial _z v_m\right ] \nonumber \\ &\quad -\, \varPi _{\mu} \left [2\epsilon ^2 \partial _z d\left (\partial _z u_a-\partial _r v_a\right )-\epsilon ^2 \partial _z v_a\right ], \end{align}

(2.6) \begin{align} p_m - \varPi _{{\mu} } p_a &= - Ca (\kappa ) + 2\epsilon ^2\left (\partial _r v_m- \partial _r u_m \partial _z d \right ) -2\epsilon ^2 \varPi _{\mu} \left (\partial _r v_a- \partial _r u_a \partial _z d \right ), \\[9pt] \nonumber \end{align}

where $Ca = \gamma /{\mu} _m U$ is the capillary number and the mean curvature $\kappa$ to $\mathcal{O}(\epsilon ^2)$ is given by

(2.7) \begin{equation} \kappa = \frac {1}{d}-\frac {\epsilon ^2\left (\partial _z d\right )^2}{2d}-\epsilon ^2 \partial _{zz}. \end{equation}

This nonlinear approximation of the full curvature is essential for capturing the onset of liquid-bridge formation (Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015), which limits the volume of mucus that can be occupied in an open airway; we will show later that nonlinear curvature also has an important effect on the extent of mucus-depleted zones at the wall. In fact, the primary reason we retain terms up to $\mathcal{O}(\epsilon ^2)$ is to consistently include these effects of nonlinear curvature in the model.

The evolution of the interface is governed by the kinematic boundary condition:

(2.8) \begin{equation} \partial _t d = v_m - u_m\partial _z d. \end{equation}

At the centreline, we have the symmetry condition,

(2.9) \begin{align} v_a = 0, \quad \partial _r u_a = 0 , \quad {\textrm{at}}\quad r = 0, \end{align}

while at the wall, we apply the non-penetration condition and impose a tangential velocity to account for mucociliary transport:

(2.10) \begin{align} v_m = 0, \quad u_m = u_c(z,t), \quad {\textrm{at}}\quad r = 1. \end{align}

The function $u_c(z,t)$ is a travelling wave, chosen to mimic the asymmetric metachronal wave of the cilia carpet:

(2.11) \begin{equation} u_c(z,t) = \begin{cases} \text{$a_f\langle \, u_c \rangle \,\sin\left (2\pi (z+ct)/\lambda _c\right )$}, \quad n\lambda _c \leq z+ct \lt n\lambda _c+\lambda _c/2,\\ \text{$a_r\langle \, u_c \rangle \,\sin\left (2\pi (z+ct)/\lambda _c\right )$}, \quad n\lambda _c+\lambda _c/2 \leq z+ct \lt n\lambda _c+\lambda _c,\\ \end{cases} \end{equation}

for $n = 0,1,2\ldots\,$ . The celerity of the wave is given by $c = f_c \lambda _c$ , where $\lambda _c$ is the metachronal wave length and $f_c$ is the cilia-beating frequency. The mean cilia velocity, $\langle \, u_c \rangle = ({1}/{\lambda _c})\int _{0}^{\lambda _c} u_c(z',t)\,{\textrm{d}}z'= ({1}/{f_c^{-1}})\int _{0}^{f_c^{-1}} u_c(z,t')\,{\textrm{d}}t'$ , is a non-zero constant because the effective forward stroke is stronger than the reverse stroke, i.e. the ratio of the constants $a_f/a_r$ is greater than unity. Note that the metachronal wave is antiplectic and propagates in the direction opposite to that of the effective cilia stroke.

This cilia boundary condition is similar to that used by Vasquez et al. (Reference Vasquez, Jin, Palmer, Hill and Forest2016) to model cilia-driven large-scale swirling flows, observed in an in vitro experiment. However, while Vasquez et al. (Reference Vasquez, Jin, Palmer, Hill and Forest2016) use a spatially invariant boundary velocity, we have used a travelling wave to better represent the metachronal cilia motion. Our boundary prescription is also closely related to that developed and validated by Bottier et al. (Reference Bottier, Peña Fernández, Pelle, Isabey, Louis, Grotberg and Filoche2017b ) and used recently by Choudhury et al. (Reference Choudhury, Filoche, Ribe, Grenier and Dietze2023) to understand the influence of viscoelasticity on mucociliary transport; this is a Navier-slip condition that contains a slip length, which when set to zero yields a Dirichlet boundary condition just like ours. We proceed with a zero-slip condition for simplicity. As demonstrated later, the cilia-induced velocity is much too slow to alter the shape of the film or the streamlines of the airflow, and so our results regarding particle deposition will not change with the addition of slip or other modifications to the precise form of the cilia boundary condition.

To ease computations, we shall later explore whether the travelling-wave cilia velocity can be replaced by just its mean value $\langle \, u_c \rangle$ ; the resulting Couette boundary condition is

(2.12) \begin{equation} v_m = 0, \quad u_m = \langle \, u_c \rangle , \quad {\textit{at}}\quad r = 1. \end{equation}

In the axial direction, we consider periodic boundary conditions. Thus, the interface profile $d(z,t)$ evolves without changing the total volume of mucus, which remains at its initial value of $(1-d_0^2) L R^2$ , where $d_0$ is the dimensionless position of the initially flat film ( $d(z,0) = d_0$ ) and $L$ is the length of the domain (see figure 1).

Equations (2.1), (2.2), (2.4)–(2.9), along with either (2.10) and (2.11), or just (2.12), form a closed system of equations. Most of the parameters are set to physiologically relevant values, listed in table 2. The initial film thickness, $1-d_0$ , is varied to study the effect of the mucus volume fraction.

Table 2.

Values of the parameters in the flow model, corresponding to the results in the main text. Some key figures are also produced for a second set of mucus-air properties in the supplementary material are available at https://doi.org10.1017/jfm.2025.10606.

2.2. Weighted-residual approach

We now average across the radial direction, using the WRIBL method, to obtain evolution equations for the interface profile $d(z,t)$ , and the flow rates of air and mucus, $Q_i(z,t)$ . The derivation is outlined, in brief, in this section; for details, the reader is referred to Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015), whom we follow closely. The implementation of this derivation using computer algebra is discussed by Hazra & Picardo (Reference Hazra and Picardo2025a ).

2.2.1. Velocity decomposition

The velocity field is decomposed into a leading-order contribution $\hat {u}_i$ and a correction $u^{\prime}_i$ that is $\mathcal{O}(\epsilon )$ . The leading term accounts for the local balance between the $\mathcal{O}(1)$ radial viscous diffusion and axial pressure-gradients, while yielding the exact flow rate $2 \pi Q_i$ . Thus, $\hat {u}_i$ is a self-similar radial profile that is parametrised by the cilia velocity, the flow rate and the film thickness:

(2.13) \begin{equation} u_i(r,z,t) = \hat {u}_i(r;u_c,d,Q_i) + u^{\prime}_i(r,z,t), \end{equation}

with

(2.14) \begin{equation} \quad \frac {1}{r}\partial _r \left (r\partial _r \hat {u}_m\right ) = A_m, \quad \frac {1}{r}\partial _r \left (r\partial _r \hat {u}_a\right ) = A_a, \end{equation}

and

(2.15) \begin{align} \hat {u}_{m}=u_{c}(z,t), \quad \hat {v}_{m} = 0\;\;\; \quad {\textit{at}}\quad r &= 1, \end{align}

(2.16) \begin{align} \partial _r \hat {u}_a = 0, \quad \hat {v}_{a} = 0\;\;\; \quad {\textit{at}}\quad r &= 0, \end{align}

(2.17) \begin{align} \partial _r \hat {u}_m = \varPi _{\mu} \partial _r \hat {u}_a, \quad \hat {u}_a = \hat {u}_m \quad {\textit{at}}\quad r &= d, \end{align}

where $A_m$ and $A_a$ are determined by the flow rates $2\pi Q_i$ :

(2.18) \begin{align} \quad \int _{0}^{d} \hat {u}_a \,r\,{\textrm{d}}r = Q_a , \quad \int _{d}^{1} \hat{u}_{m}r{\textrm{d}}r = Q_m. \end{align}

Once $d$ and $Q_i$ are obtained by solving the WRIBL equations, given in the next sub-section, then $\hat {u}_i$ can be calculated to obtain the leading axial velocity profile at any position in the airway. Furthermore, the leading radial-velocity profile can be obtained by integrating the continuity equation (2.1):

(2.19) \begin{align} \quad \hat {v}_a = -\frac {1}{r}\int _{0}^{r} \partial _z\hat {u}_a r\,{\textrm{d}}r ,\quad \hat {v}_m = \frac {1}{r}\int _{r}^{1} \partial _z\hat {u}_m r\,{\textrm{d}}r. \end{align}

In this manner, the leading-order, incompressible velocity field $(\hat {u}_i,\hat {v}_i)$ is obtained. It is used to (i) evolve particles in the air phase, and (ii) visualise the flow by plotting streamlines obtained from contours of the streamfunction $\varPsi _i$ :

(2.20) \begin{align} \quad \varPsi _a = \int _{0}^{r} \hat {u}_a \,r \,{\textrm{d}}r, \quad \varPsi _m = \varPsi _a|_d+\int _{d}^{r} \hat {u}_m \,r\, {\textrm{d}}r. \end{align}

Regarding the velocity corrections $u^{\prime}_i$ , (2.18) implies that they must satisfy the following guage conditions, which will be used to eliminate $u^{\prime}_i$ from the weighted-averaged equations,

(2.21) \begin{align} \quad \int _{0}^{d} u^{\prime }_a \,r\,{\textrm{d}}r = 0 , \quad \int _{d}^{1} u^{\prime }_m \,r\,{\textrm{d}}r = 0. \end{align}

2.2.2. Weighted residuals

Denoting the boundary layer equations for the two phases in (2.2) by $\textit{BLE}_i$ , we obtain averaged equations by evaluating the residual $\langle \textit{BLE}|w \rangle$ where the inner product is defined as $\langle p|q \rangle = \varPi _{\mu} \int _0^d{p_a q_a r\,{\textrm{d}}r} + \int _d^1{p_m q_m r\,{\textrm{d}}r}$ and $w_i$ are weight functions. If naive radial averaging is performed with $w_i=1$ , then the velocity correction $u^{\prime}_i$ will arise in the residual via the leading-order, radial viscous diffusion term. One will then have to either calculate $u^{\prime}_i$ or suffer prominent errors that can produce unphysical solutions (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). The ingenuity of the WRIBL approach (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012; Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2013) is to use a weight function which eliminates the correction from the residual (to $\mathcal{O}(\epsilon ^2))$ .

Multiple choices for the weights are possible. We follow Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2013, Reference Dietze and Ruyer-Quil2015) and use two combinations. First, we choose $w_i$ so that the pressure-gradient term is cancelled out and an evolution equation for $Q_i$ is obtained. The defining equations for this weight function are

(2.22) \begin{align} \quad \frac {1}{r}\partial _r \left (r\partial _r {w}_m\right ) = 1, \quad \frac {1}{r}\partial _r \left (r\partial _r {w}_a\right ) = C_a, \end{align}

with the homogeneous version of the $\hat {u}_i$ boundary conditions:

(2.23) \begin{align} {w}_{m}|_{r=1}=0, \!\quad \partial _r {w}_a|_{r=0} = 0, \!\quad \partial _r {w}_m|_{r=d} = \varPi _{\mu} \partial _r {w}_a|_{r=d}, \!\quad {w}_a|_{r=d} = {w}_m|_{r=d}. \end{align}

Here, $C_a$ in (2.22) follows from applying

(2.24) \begin{align} \int _{0}^{d} w_a r\,{\textrm{d}}r = -\int _{d}^{1} w_m r\,{\textrm{d}}r. \end{align}

The second weight function, denoted by $\tilde {w}_i$ , is chosen to yield a diagnostic equation for the pressure at the interface, $p_i|_d$ . This equation will be used to enforce axial pressure boundary conditions. The defining equations for $\tilde {w}_i$ are

(2.25) \begin{align} \quad \frac {1}{r}\partial _r \left (r\partial _r \tilde {w}_m\right ) = 1, \quad \frac {1}{r}\partial _r \left (r\partial _r \tilde {w}_a\right ) = \tilde {C}_a, \end{align}

with the same boundary conditions as $w_i$ , but with

(2.26) \begin{align} \int _{0}^{d} \tilde {w}_a r\,{\textrm{d}}r = \int _{d}^{1} \tilde {w}_m r\,{\textrm{d}}r, \end{align}

instead of (2.24).

On applying the velocity decomposition (2.13), along with (2.21), the corrections $u^{\prime}_i$ cancel out from the residual $\langle \textit{BLE}|w \rangle$ at $\mathcal{O}(\epsilon ^2)$ , provided we neglect $\mathcal{O}({\textit{Re}}_i \epsilon ^2)$ inertial terms, as is appropriate for low to moderate Reynolds numbers (Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015).

Note that our choice of weight functions differs from that of Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015). Our final WRIBL equations, listed later, have the same form as those of Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015), and the coefficients also match after rescaling to account for the difference in weight functions. We have also checked that the predictions of our WRIBL equations match those of Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015), as illustrated by figure 15 in the Appendix.

2.2.3. WRIBL equations

Two exact mass-conservation equations follow from integrating the continuity equation (2.1) in the radial direction, in each phase, and using the kinematic boundary condition (2.8):

(2.27) \begin{align} d \partial _t d &= \partial _z Q_m, \end{align}

(2.28) \begin{align} \partial _z Q_m &+ \partial _z Q_a = 0, \end{align}

where, alternately, the latter could be replaced by $ d\partial _t d = -\partial _z Q_a$ .

The evolution equation for the flow rate follows from evaluating the weighted residual with respect to $w_i$ :

(2.29) \begin{align} \frac {{\mu} _i}{{\mu} _m}{\textit{Re}}_i & \left (S_{\textit{ij}}\partial _t Q_j+S_{\textit{ic}}\partial _t u_c+F_{\textit{ijk}} Q_j{\partial _z Q_k} +F_{ijc} Q_j\partial _z u_c+F_{\textit{icj}} u_c\partial _z Q_j+ F_{\textit{icc}}u_c\partial _zu_c \right. \nonumber \\ &\quad \left. +\, G_{\textit{ijk}}Q_jQ_k\partial _z{d}+G_{\textit{icj}}u_cQ_j\partial _z{d}+G_{\textit{icc}}u_cu_c\partial _z{d} \right)\nonumber \\ &=-Ca\left (\partial _z\kappa +\partial _z\kappa _p\right )Iw_a+\varPi _{{\mu} }C_aQ_a+Q_m+J_{j}Q_{j}(\partial _z {d})^2+J_{c}u_{c}(\partial _z {d})^2\nonumber \\ &\quad +\, K_{j}\partial _z Q_{j}\partial _z {d}+K_{c}\partial _z u_{c}\partial _z {d} +L_{j}Q_{j}\partial _z^2 {d}+L_{c}u_{c}\partial _z^2 {d}\nonumber \\&\quad +\, M_j\partial _z^2Q_j+M_c\partial _z^2u_c-u_c\partial _r{w_m}|_{1}+\varPi _{{\mu} } G_{a}\sin(\omega t)Iw_a, \end{align}

where Einstein summation notation has been used on the phase index. The coefficients, $S_{\textit{ij}},S_{\textit{ic}},F_{\textit{ijk}},\ldots ,M_c$ , are functions of $d$ alone. The integrals of $\hat {u}$ and $w_i$ , which yield these coefficients, are involved and the corresponding calculations are performed using the Python symbolic-computing library SymPy (Meurer et al. Reference Meurer2017).

Note that, for convenience of notation, we have written the WRIBL (2.29) after rescaling the axial coordinate with $R$ ; thus, the length-scale ratio $\epsilon$ does not appear here. This new axial scaling will be used to present all subsequent equations and figures.

Equations (2.27)–(2.29) form a closed system for $d$ , $Q_m$ and $Q_a$ . Integrating (2.28) shows that the total flow rate, $Q_t(t) = Q_a+Q_m$ , is spatially invariant. Therefore, given $Q_t$ as an input, one can replace $Q_a$ by $Q_t-Q_m$ and solve (2.27) and (2.29) for $d$ and $Q_m$ . Here, however, we do not impose $Q_t$ . Rather, we impose the background pressure-gradient, so that the net pressure difference across the computational domain, of length $L$ , must be $G_{a}L\,\sin(\omega t)$ in the air and zero in the mucus. The deviation from this applied pressure-gradient, $p_i$ , must therefore integrate to zero over the domain:

(2.30) \begin{equation} \int _{0}^{L/R } \partial _z p_i\, {\textrm{d}}z = 0. \end{equation}

To apply this pressure constraint, we require the diagnostic equation for pressure, which is obtained by evaluating the weighted residue with respect to $\tilde {w}_i$ :

(2.31) \begin{align} \frac {{\mu} _i}{{\mu} _m}{\textit{Re}}_i & \left(\tilde {S}_{\textit{ij}}\partial _t Q_j+\tilde {S}_{\textit{ic}}\partial _t u_c+\tilde {F}_{\textit{ijk}} Q_j\partial _z Q_k +\tilde {F}_{ijc} Q_j\partial _z u_c+\tilde {F}_{\textit{icj}} u_c\partial _z Q_j+\tilde {F}_{\textit{icc}}u_c\partial _zu_c \right. \nonumber \\ &\quad\quad\quad\,\,\,\, \left. +\,\tilde {G}_{\textit{ijk}}Q_jQ_k\partial _z{d}+\tilde {G}_{\textit{icj}}u_cQ_j\partial _z{d}+\tilde {G}_{\textit{icc}}u_cu_c\partial _z{d}\right)\nonumber \\ & =-2\varPi _{{\mu} }\partial _z p_a |_d\tilde {I}w_a+Ca\left (\partial _z\kappa +\partial _z\kappa _p\right )\tilde {I}w_a+\varPi _{{\mu} }\tilde {C}_aQ_a+Q_m+\tilde {J}_{j}Q_{j}(\partial _z {d})^2\nonumber \\ &\quad +\, \tilde {J}_{c}u_{c}(\partial _z {d})^2 + \tilde {K}_{j}\partial _z Q_{j}\partial _z {d}+\tilde {K}_{c}\partial _z u_{c}\partial _z {d} +\tilde {L}_{j}Q_{j}\partial _z^2 {d}+\tilde {L}_{c}u_{c}\partial _z^2 {d}\nonumber \\ &\quad +\, \tilde {M}_j\partial _z^2Q_j+\tilde {M}_c\partial _z^2u_c-u_c\partial _r{\tilde {w}_m}|_{1}+\varPi _{{\mu} } G_{a}\sin(\omega t)\tilde {I}w_a. \end{align}

As is customary when using thin-film models to describe droplets, rivulets and other situations involving regions devoid of liquid (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Ghatak, Khanna & Sharma Reference Ghatak, Khanna and Sharma1999; Ruyer-Quil et al. Reference Ruyer-Quil, Bresch, Gisclon, Richard, Kessar and Cellier2023), we employ a precursor-film approach to capture the mucus-depleted zones of the airway. To wit, we introduce a disjoining pressure via the term $\kappa _p = {h^6_0}/{(1-d)^9}$ in (2.29) and (2.31). With $h_0 = 5 \times 10^{-4}$ , we obtain a precursor film thickness below 0.01 ( $d\gt 0.99$ ) in the depleted zones; the fluid drained out of these zones accumulates in annular humps whose shape, which we have verified, is very close to that of equilibrium unduloids (of the same volume). Reducing the precursor film’s thickness even further does not affect our results but significantly increases the computational cost.

The coefficients of (2.29) and (2.31) are stored in text files and imported when performing numerical simulations, as done in the Jupyter notebook associated with figure 15 in the Appendix.

2.2.4. One-way coupled WRIBL equations

Given that air has a much smaller viscosity and density than mucus, it is natural to check whether the simplified model obtained in the limit of $\varPi _{\rho }$ , $\varPi _{{\mu} } \ll 1$ provides a good approximation of the two-phase flow. In this limit, the mucus film evolves independently of the air, with a stress-free surface. Such a passive-core approximation is routinely made when studying gas–liquid flows in which the gas flow is weak. While this approximation is certainly inappropriate for the upper airways, where the turbulent airflow is strong enough to rip off mucus droplets from the film, it could work well in the middle airways where the airflow is much weaker.

The WRIBL model for the mucus phase in the passive-core limit was analysed by Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015). Here, we additionally derive the WRIBL equation for the air phase, because we need the air velocity field to study aerosol transport. In the limit of infinite viscosity contrast, the air will experience the mucus film as a moving solid boundary. Thus, while the mucus film evolves independently, the airflow is modulated by the deformation of the film, i.e. we have a one-way coupled flow. Because the mucus velocity is typically much slower than that of the air, the major effect of the film is to alter the conduit for airflow.

The decoupled mucus flow is governed by the domain equations (2.1) and (2.2), with $i = m$ , the wall boundary-condition (2.10), the kinematic condition (2.8), and the free surface conditions obtained by setting $\varPi _{\mu} = 0$ in (2.5) and (2.6). (The smallness of $\varPi _{\rho }$ is used implicitly when we set $\varPi _{\mu} p_a$ to zero, in (2.6); $\varPi _{\rho }$ controls the relative strength of air pressure variations, in cases where inertia is significant.) We then apply the WRIBL method to the mucus phase by performing the weighted integral from $r=d$ to $r = 1$ . The leading-order mucus velocity is determined by the equations for $\hat {u}_m$ in (2.14)–(2.18), but now with $\varPi _{\mu} = 0$ . Similarly, the weight function is defined by the equations for $w_m$ in (2.22) and (2.23), with $\varPi _{\mu} = 0$ . The resulting WRIBL model for the mucus film is as follows:

(2.32) \begin{align}{\textrm{d}}\partial _t d &= \partial _z Q_m,\end{align}

(2.33) \begin{align} {\textit{Re}}_m & \left(\overline {S}_{\textit{mm}}{\partial _t Q_m}+\overline {S}_{mc}{\partial _t u_c}+\overline {F}_{mmm} Q_m\partial _z Q_m +\overline {F}_{\textit{mmc}}Q_m\partial _z u_c+\overline {F}_{\textit{mcc}}u_c\partial _z u_c \right. \nonumber \\ &\quad \left. +\, \overline {F}_{mcm} u_c\partial _z Q_m+\overline {G}_{mmm}Q_mQ_m\partial _z{d}+\overline {G}_{mcm}Q_m u_c\partial _z{d}+\overline {G}_{\textit{mcc}}u_c u_c\partial _z{d}\right)\nonumber \\&\quad =Ca\left (\partial _z\kappa +\partial _z\kappa _p\right )\overline {Iw}_m+Q_m+\overline {J}_{m}Q_{m}(\partial _z {d})^2+\overline {J}_{c}u_{c}(\partial _z {d})^2+\overline {K}_{m}\partial _z Q_{m}\partial _z {d}\nonumber \\ &\quad +\, \overline {K}_{c}{\partial _z u_{c}}{\partial _z {d}} +\overline {L}_{m}Q_{m}{\partial _z^2 {d}}+\overline {L}_{w}u_{c}{\partial _z^2 {d}}+\overline {M}_m{\partial _z^2Q_m}+\overline {M}_c{\partial _z^2u_c}-u_c\partial _r{\overline {w_m}}|_{1}. \end{align}

Equations (2.32) and (2.33) must be solved simultaneously to determine $Q_m$ and $d$ , which will then be used as inputs to calculate $Q_a$ from the WRIBL equation in the air phase. The leading-order air velocity $\hat {u}_a$ is determined by (2.14) and (2.18) along with the boundary conditions $\partial _r \hat {u}_a|_0 = 0$ and $\hat {u}_a|_d = \hat {u}_m|_d$ . The weight function $w_a$ satisfies the homogeneous version of these boundary conditions, in addition to (2.22). Performing the weighted integral over the air phase, from $r=0$ to $r = d$ , yields the WRIBL equation for the airflow:

(2.34) \begin{align} {\textit{Re}}_a & \left(\overline {S}_{\textit{aa}}{\partial _t Q_a}+\overline {S}_{\textit{ad}}{\partial _t u_d}+\overline {F}_{aaa} Q_a{\partial _z Q_a} +\overline {F}_{aad}Q_a{\partial _z u_d} + \overline {F}_{ada} u_d{\partial _z Q_a} \right. \nonumber \\ & \quad \left. +\, \overline {F}_{\textit{add}}u_d\partial _z u_d+\overline {G}_{aaa}Q_aQ_a{\partial _z{d}}+\overline {G}_{ada}u_dQ_a{\partial _z{d}}+\overline {G}_{\textit{add}}u_du_d{\partial _z{d}}\right)\nonumber \\ &=-\partial _z p_a|_d\overline {Iw}_a+\overline {C}_aQ_a+\overline {J}_{a}Q_{a}(\partial _z {d})^2+\overline {J}_{d}u_{d}(\partial _z {d})^2+\overline {K}_{a}{\partial _z Q_{a}}{\partial _z {d}}+\overline {K}_{d}{\partial _z u_{d}}{\partial _z {d}}\nonumber \\ &\quad +\,\overline {L}_{a}Q_{a}{\partial _z^2 {d}}+\overline {L}_{d}u_{d}{\partial _z^2 {d}}+\overline {M}_a{\partial _z^2Q_a}+\overline {M}_d{\partial _z^2u_d}+G_{a}sin(\omega t)\overline {Iw}_a-u_d\partial _r{\overline {w_a}}|_{d}. \end{align}

This equation, which has the air flow rate $Q_a$ and the pressure at the interface $p_a|_d$ as unknowns, must be solved along with the overall mass balance equation (2.28).

2.2.5. Numerical solution

To solve the fully coupled WRIBL equations, we first use (2.28) to substitute $Q_m = Q_t(t)-Q_a$ in (2.29) and (2.31). The latter pressure equation is then integrated over the domain and the pressure constraint (2.30) is applied to obtain an ordinary differential equation (ODE) for $Q_t(t)$ of the form

(2.35) \begin{equation} {d}_t Q_t = \varPhi (d,Q_a,Q_t,u_c). \end{equation}

Equations (2.27), (2.29) and (2.35) thus become a closed system for $d$ , $Q_a$ and $Q_t$ .

Numerical simulations are performed by discretising space using a second-order central-difference scheme. The resulting system of ODEs is integrated in time using the stiff, adaptive time-stepping, LSODA solver (Hindmarsh & Petzold Reference Hindmarsh and Petzold1995) provided by the function solve_ivp, of the Python library SciPy (Virtanen et al. Reference Virtanen2020).

Turning to the one-way coupled WRIBL model, the equations for the mucus film, namely (2.32) and (2.33), are solved first to obtain $d$ and $Q_m$ . Here too, we employ second-order central-differencing in space, followed by the LSODA method for stepping in time. To obtain $Q_a$ , we need only solve for $Q_t$ , since $Q_a = Q_t(t)-Q_m$ ; on substituting this relation into the one-way coupled WRIBL equation for air, (2.34), integrating over the domain and applying the pressure constraint (2.30), we obtain an ODE for $Q_t$ . So, rather than simultaneously solving two partial differential equations (PDEs) and one ODE, as required by the fully coupled WRIBL model, we solve just two PDEs first and then separately solve a single ODE.

All our simulations are performed with periodic boundary conditions in the axial $z$ direction. The length of the computational domain $L$ is chosen to match the wavelength of the fastest growing mode of the Rayleigh–Plateau instability, $\varLambda _{\textit{RP}} = 2\pi 2^{1/2} d_0 R$ . (This inviscid prediction of Rayleigh (Reference Rayleigh1892) works very well even for viscous films (Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015)). An evenly spaced spatial grid of 500 points is found to be sufficient for obtaining grid-independent solutions.

2.3. Particle motion

To track the motion of airborne particles, we adopt the point-particle approximation and use the Maxey–Riley equation, simplified for tiny dense particles (Ravichandran et al. Reference Ravichandran, Deepu and Govindarajan2017). In addition, we include Brownian forces to describe a wide range of particles, from Brownian to inertial. The particle’s position $\boldsymbol{X}_p$ and velocity $\boldsymbol{V}_p$ are thus governed by the following stochastic differential equations (SDEs), which are often used for studying particle transport in the lungs (Tsuda et al. Reference Tsuda, Henry and Butler2013):

(2.36) \begin{align} {\textrm{d}}\boldsymbol{X}_p &= \boldsymbol{V}_p\, {\textrm{d}}t, \end{align}

(2.37) \begin{align} St\,{\textrm{d}}\boldsymbol{V}_p &= \left [\boldsymbol{u}(\boldsymbol{X}_p,t) {-\boldsymbol{V}_p}\right ]{\textrm{d}}t+\sqrt {\frac {2}{Pe}}\, {\textrm{d}}\boldsymbol{W}, \end{align}

where $\boldsymbol{u}$ is the air-velocity vector which must be obtained at the instantaneous position of the particle, $\boldsymbol{X}_p$ . Equation (2.37) equates the rate of change in the particle’s momentum to the sum of the Stokes drag force (proportional to the slip velocity, $\boldsymbol{u}-\boldsymbol{V}_p$ ) and the Brownian force. The latter is modelled by ${\textrm{d}}\boldsymbol{W}$ , a vector of independent increments of the Wiener process (Gardiner Reference Gardiner2009). The impact of Brownian forces on the particle’s motion varies inversely with the Péclet number ( $\textit{Pe}$ ), which is the ratio of the time scales of Brownian diffusion to convection by the airflow. The importance of inertial effects increases with the Stokes number ( $St$ ), which is the ratio of the time scale of the particle’s inertial relaxation to the time scale over which the flow changes (as viewed by the convected particle). In our study, $\textit{Pe}$ and $St$ are not independent but vary together with the particle diameter $d_p$ , as follows:

(2.38) \begin{equation} Pe = \frac {R U^\star }{D} \; \textrm{with} \; D = \frac {k_B T}{3\pi {\mu} _a d_p} \ \textrm{and}\ St = \frac {\frac {\pi }{6} \rho _p d_p^3}{\zeta } \frac {U^\star }{R}= \frac {d_p^2 \rho _p}{18 {\mu} _a} \frac {U^\star }{R}, \end{equation}

where $D$ is the Einstein diffusivity of the particle, $k_B$ is the Boltzmann constant, $T = 298 \,\textrm{K}$ is the temperature, $\rho _p$ is the particle’s density and $\zeta = 3 \pi {\mu} _a d_p$ is the Stokes friction coefficient. The air velocity scale $U^\star$ is chosen to be $U_{{max}}/2$ , where $U_{{max}} = {max}_t [2({ Q_a|_{d_{{min}}}})/{d_{min}^2} ] U$ is the maximum, over time, of the cross-section averaged air velocity below the mucus hump (the $z$ location where $d = d_{{min}}$ ); it yields an estimate of the smallest convective time scale during respiration. We fix the particle’s density to that of water ( $\rho _p = 1000$ $\textrm{kg m}^{-3}$ ) and vary the diameter from $10^{-3}$ to 50 $\unicode{x03BC}\textrm{m}$ . Both $\textit{Pe}$ and $St$ increase strongly with $d_p$ (see table 3) so that the motion of the smallest particles is nearly Brownian, while that of the largest is influenced by inertial effects, such as migration across curved streamlines. Intermediate-sized particles are affected the least by these forces and tend to follow the flow like tracers.

Table 3.

Range of particle sizes considered in this study and the corresponding time scales and non-dimensional numbers.

In solving (2.36)–(2.37), we use a two-step exponential integration scheme, which enables accurate simulations of small $St$ particles (Ireland et al. Reference Ireland, Vaithianathan, Sukheswalla, Ray and Collins2013). The Wiener increment ${\textrm{d}}\boldsymbol{W}$ is updated only once every time step (it stays the same for both half-steps); thus, the Brownian term is integrated using the Euler–Marujama scheme (Öttinger Reference Öttinger1996). This strategy, of using higher-order time-stepping for the drift terms of an SDE, is helpful when the relative importance of Brownian forces varies widely; it has been used, for example, to study Brownian elastic chains in turbulent flows (Bagheri et al. Reference Bagheri, Mitra, Perlekar and Brandt2012; Singh et al. Reference Singh, Picardo and Ray2022).

The solution of the WRIBL model yields the air velocity field in axisymmetric cylindrical coordinates. We transform this field to three-dimensional Cartesian coordinates and then evolve particles in three dimensions. The air velocity at the particle location is obtained from its value on the computational grid via linear interpolation. We consider the dilute limit and neglect inter-particle collisions. Thus, the motion of each particle is independent of all others in the airway; we use this fact to parallelise the numerical integration of the particles to accelerate simulations. We track up to 2000 particles in a given flow, for up to 30 breathing cycles.

Adopting the simplest approach to studying deposition, we assume that a particle deposits when its surface makes contact with the film (Guha Reference Guha2008). We therefore monitor the radial distance (along the $r$ coordinate of the airway) between the centre of mass of the particle and the mucus–air interface; if this distance falls below $d_p/2$ , then we assume that the particle has collided with the mucus film and become instantly entrapped; the trapped particle is no longer evolved. Strictly speaking, we should monitor the shortest distance between the particle and the interface, but this would be computationally expensive. So, noting that axial variations of the film occur over distances much larger than the particle’s radius, we replace the shortest distance by the radial distance. In keeping with our use of a precursor film to represent mucus-depleted sections of the wall, we treat all particles that collide with the precursor film (of thickness less than $0.01R$ ) as having deposited on the wall.

This model of particle transport, while simple, captures the key physical aspects of deposition in a mucus-lined airway: (i) size-dependent Brownian and inertial forces; (ii) possibility of deposition on mucus-depleted zones of the wall; and (iii) influence of spatially non-uniform airflow, due to the presence of protruding mucus humps. There are of course various directions in which the particle model can be improved. In particular, the largest particles considered here are 1/8 times the diameter of the airway and are thus not well suited to be treated as point particles. Nonetheless, we include these particles in our study, with the aim of providing a first approximation to their deposition behaviour. Indeed, the influence of finite-size effects (such as shear-gradient-induced lift and wall-induced lift) on the particle’s motion scale with the fluid to particle density ratio (Guha Reference Guha2008; Liu et al. Reference Liu, Xue, Sun and Hu2016); here, $\rho _p/\rho _a = 10^{-3}$ and so the effects of finite size will be weak in comparison to the effects of particle inertia (such as centrifugation). The latter will be shown to play an important role in the presence of a non-uniform mucus film.

3. How cilia and air affect the mucus film

We shall first examine how the mucus film – its distribution and flow – is affected by cilia and air. Then, in § 4, we will quantify the extent of mucus-depleted zones. The understanding of the flow gained in these sections will aid our study of particle deposition, which will be taken up in § 5.

All our simulations use the parameters given in table 2, which are typical for mucus and air in a middle airway (we focus on generations 14–16). Our results are not sensitive to the precise values of the parameters, as demonstrated in the supplementary material which reproduces select figures for a different value of the viscosity ratio and interfacial tension.

Throughout this work, we set the length of the computational domain to match the wavelength of the fastest growing mode of the Rayleigh–Plateau instability, i.e. $L = \varLambda _{\textit{RP}}= 2\pi 2^{1/2} d_0 R$ (Rayleigh Reference Rayleigh1892; Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015). Furthermore, we only consider films with an initial thickness $1-d_0$ sufficiently small to remain open without forming a liquid bridge, i.e. with $(1-d_0^2) \varLambda _{\textit{RP}}/R \lt 1.73$ (Everett & Haynes Reference Everett and Haynes1972; Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015). In the present section, representative results for the air–mucus flow are shown for the value of $d_0 = 0.885$ . The simulations are initialised with a slightly perturbed interface, $d = d_0 + 10^{-3}\,\textrm{cos}(2 \pi z R/\varLambda _{\textit{RP}})$ . Random perturbations will be applied later, in § 4, while analysing the extent of mucus-depleted zones.

3.1. Cilia translates the mucus film

The effect of cilia on the film profile is illustrated in figure 2. Snapshots of the evolution of the interface with ciliary forcing (via the metachronal velocity $u_c$ imposed at the wall, as per (2.11)) are compared with those without ciliary forcing (setting $u_c = 0$ ), in panels (a–c) and (d–f). The same unduloid-shaped hump is formed in both cases, but the hump translates in the presence of cilia. The growth of the hump, visualised by the time trace of the minimum position of the interface in figure 2(g), is indeed nearly identical in the two cases. In this and all subsequent plots, we present the temporal evolution in multiples of the non-dimensional breathing time-period $T_b$ (where the dimensional time period, $T_b R/U$ is 1 s).

The cilia-induced translation shows up clearly in figure 2(h), which traces the location of the centre-of-mass of the film, $X_{cm}$ ; the translation is near-linear with a constant speed equal to the mean cilia speed $\langle u_c \rangle$ . When this linear translation is subtracted out, as done in the inset, we find that tiny oscillations are superimposed on the linear translation; the oscillations are caused by respiratory airflow as evidenced by their presence in the case without cilia (for which we use $\langle u_c \rangle = 0$ in the inset).

These air-induced oscillations are also seen in the root-mean-square averaged velocity of mucus, $\langle u_m\rangle _{{rms}} = [{L^{-1}}\int _0^L{ ({2 Q_m}/({1-d^2}) )^2 \,{\textrm{d}}z} ]^{1/2}$ , which is plotted in figure 2(i). Here, the prominent peak in $\langle u_m\rangle _{{rms}}$ corresponds to the relatively rapid capillary-driven flow of mucus that produces the hump. After the hump has formed (and the capillary forces due to axial and radial curvatures have balanced out), the interface profile remains unchanged while the cilia slowly transports the film.

Figure 2.

Cilia translates the mucus film without altering the capillary-driven dynamics and the consequent emergence of humps and depleted zones. The top two rows present snapshots of the evolution of the film, (a–c) with cilia and (d–f) without cilia, i.e. with $u_c = 0$ . Panel (g) compares the growth of the hump by tracing the evolution of the minimum of $d(z,t)$ . Panel (h) compares the time-traces of the centre of mass of the entire film, clearly showing that cilia causes the film to translate; the inset subtracts out the net translation $\langle u_c \rangle t$ and reveals tiny lateral air-induced oscillations. Panel (i) presents the time-trace of the spatially averaged root-mean-square mucus velocity, normalised with the mean cilia velocity; the inset zooms into the curve for the no-cilia case.

Figure 3.

Metrachronal cilia velocity (2.11) has the same effect on the flow as the constant cilia velocity (2.12), as evidenced by (a, c) the kymographs of the evolution of the film thickness $1-d(z,t)$ and (b, d) the snapshots of the streamlines in the air and mucus. The snapshot of the streamlines corresponds to the time of the second line profile drawn in the kymographs ( $t \approx 15.6 \,T_b$ ).

Figure 4.

Comparison of the time traces of (a) the minimum position of the interface, (b) the centre-of-mass of the film, and (c) the root-mean-square mucus velocity, for the metrachronal and constant cilia velocities (given by (2.11) and (2.12), respectively). The curves overlap in each panel, showing perfect agreement.

Since hump formation is a one-time event, we will focus on the long term motion of the film for our study of particle deposition. The space–time kymograph of the interface profile in figure 3(a) emphasises that this long-term motion is dominated by simple translation. In fact, if we replace the asymmetric metachronal wave form of $u_c$ in (2.11) by just its mean, as given by the constant velocity Couette condition in (2.12), we obtain a near-identical translatory motion of the interface, as shown in figure 3(c). The streamlines in the air and mucus phases are also found to match very well (compare figures 3 b and 3 d). Further evidence for this equivalence is provided in figure 4 in which the time-traces of $d_{min}$ , $X_{cm}$ , and $\langle u_m\rangle _{{rms}}$ for the constant cilia velocity are seen to match those for the metrachronal cilia velocity. Thus, the metachronal waveform of $u_c$ plays no role in the evolution of the interface and the net mucociliary transport – only the mean $\langle u_c \rangle$ matters (we have checked that this result is independent of spatial resolution by increasing the grid points from 300 to 1024). Hence, in all subsequent results, we shall use the constant cilia velocity boundary condition (2.12), both for simplicity and ease of computation (without the relatively rapid oscillations in $u_c$ , we can take much larger time steps while integrating the WRIBL equations).

Of course, our conclusion regarding the unimportance of the metachronal waveform is valid only for the slow evolution of the film over long times scales; indeed, this result is unsurprising given that the cilia-beating time period is approximately three orders of magnitude smaller than the time scale with which mucus is transported across the airway (table 1). Another caveat is that we have considered a Newtonian mucus film, in which the dominant mechanism of momentum transport from the ciliary boundary is viscous diffusion, which effectively filters out the rapid oscillations of $u_c$ . In contrast, for strongly viscoelastic mucus, present under diseased conditions, Choudhury et al. (Reference Choudhury, Filoche, Ribe, Grenier and Dietze2023) have shown that the metachronal waveform of $u_c$ produces a substantial reduction of the net mucociliary transport (in a flat, non-deforming mucus film). However, for healthy mucus, which is not as viscoelastic, this effect is marginal; indeed, the flow rate of a flat Newtonian film is shown to depend only on $\langle u_c \rangle$ (Choudhury et al. Reference Choudhury, Filoche, Ribe, Grenier and Dietze2023).

The cilia-induced translation of mucus is much slower than the respiratory airflow, as is evident in figure 2(h): the interface hardly moves in a single breath, i.e. in a time interval of $T_b$ (see also the corresponding time and velocity scales in table 1). So, the interface will appear to be stationary to air-borne particles which typically travel back and forth across the domain with each breath. Indeed, we see in figure 3(b) that the air streamlines bend around the mucus hump as if it were stationary. So, while ciliary transport of mucus is essential for the ultimate evacuation of trapped particles from the lungs, it cannot play a role in the deposition of particles on the mucus; mucociliary transport is simply too slow relative to the motion of airborne particles to produce any collisions that would not occur without cilia.

3.2. Airflow does not alter the distribution of the mucus film

Next, we consider the effect of respiratory airflow on the film. The small oscillations in figure 2(i) show that the flow of mucus is marginally modulated by the airflow. This effect would obviously not be captured by the one-way coupled WRIBL model, in which the flow of air does not affect that of mucus. However, aerosol transport and deposition depend crucially on the mucus film’s distribution (shape and size of humps and depleted zones) and not on its internal flow. Moreover, because the flow of mucus is much slower than that of air, the film will influence the airflow primarily via changes in its profile. The interface profile $d(z,t)$ , and the flow of air around it, might well be captured by the one-way model even though the flow inside the mucus film is not exactly reproduced.

Figure 5.

Comparison of the interface evolution predicted by the fully coupled and the one-way coupled WRIBL models; in the latter, air does not affect the flow mucus film. (a) Evolution of the minimum interface position (main panel) and the centre of mass (inset). Kymographs of the film thickness are compared in panels (b–c).

Figure 6.

Comparison of the evolution of averaged velocities in (a) mucus and (b) air, as predicted by the fully coupled and the one-way coupled WRIBL models. Panel (a) presents the spatially-averaged root-mean-square mucus velocity; note the absence of the small air-induced oscillations in prediction of the one-way coupled model, which ignores the influence of air on the mucus film. Panel (b) presents the radial average of the air velocity just below the hump, i.e. below the thickest part of the interface given by $d_{min}$ .

Figure 5 shows that the one-way coupled model does indeed capture the evolution of the interface almost perfectly. The time traces of $X_{cm}$ and the minimum of $d$ practically overlap (figure 5 a), and the kymographs exhibit the same translating hump (figure 5 b,c).

Is the airflow reproduced as accurately by the one-way coupled model? Figure 6 compares time traces of (a) the root-mean-square mucus velocity and (b) the radially averaged axial velocity below the hump, obtained from the fully coupled and the one-way coupled models. We see that though the mucus flow is not entirely reproduced – the small air-induced oscillations are missing – the airflow does appear to be accurately predicted by the one-way coupled model.

Figure 7 compares the predictions of the airflow is more detail. Panels (a) and (b) zoom into a few oscillations of the long-term air and mucus flow. The remaining panels present snapshots of the instantaneous streamlines, for five selected time points within one breathing cycle, $t_1 {-} t_5$ in figure 7(a). Comparing the streamline plots in the left column (fully coupled model) with those on the right (one-way coupled model), we see that the air streamlines are perfectly captured by the one-way model except for the tiny duration of time between inhalation and exhalation when the airflow is reversing (see figures 7 g and 7 h corresponding to time $t_3$ ). For most of the breathing cycle, the airflow is orders of magnitude faster than the mucus flow (compare the y-axis ranges of figures 7 a and 7 b), and so the air flows past the mucus film as if it were a stationary solid wall. However, when the air velocity is reversing, its magnitude becomes comparable to that of the mucus velocity for a short period of time; this is when the flows in the two fluids are strongly coupled. However, the corresponding time duration is extremely small, as evidenced by the snapshots at times $t_2$ and $t_4$ just before and after $t_3$ (see figure 7 a). At both these times, the air streamlines again flow around the mucus as if it were stationary, and the one-way coupled model captures the airflow very well (figures 7 e and 7 f for $t_2$ , 7 i and 7 j for $t_4$ ).

Thus, from the perspective of particle transport and deposition, the one-way coupled model is well suited to describe the flow of air: particles will hardly move during the time for which the airflow is inaccurately predicted by the one-way model, because not only is the corresponding time duration very small, but also the airflow speed is nearly zero. We will therefore use the one-way coupled WRIBL model to obtain the interface profile and the air flow field for our analysis of particle transport in § 5.

Figure 7.

Detailed comparison of the air and mucus flows predicted by the fully coupled and one-way coupled WRIBL models; the latter ignores the influence of air on the mucus film. Panels (a) and (b) are temporal zooms of figures 6(b) and 6(a), respectively, showing the temporal variation of the mean air velocity (below the hump) and the mean mucus velocity. The line-legend is the same as that in figure 6. Five time-points are chosen, labelled $t_1$ to $t_5$ , corresponding to times when the airflow is opposite to cilia transport ( $t_1$ ), is in the same direction as cilia transport ( $t_2$ ), and is reversing ( $t_3 {-} t_5$ ). The streamlines in the air and mucus at these five times are compared in panels (c–l) for the fully coupled (left column, with inset time labels) and one-way coupled (right column) models.

4. Extent of mucus-depleted zones on the wall

We have seen, in the previous section, that the mucus film ultimately forms an unduloid-like hump that spans a portion of the wall, leaving the rest of the wall without mucus and exposed to particles (figure 2 a–c). We now examine how the extent of these mucus-depleted zones varies on increasing the volume or, equivalently, the initial thickness $1-d_0$ of the film.

We have performed numerous simulations of randomly perturbed films, with different initial thicknesses ( $1-d_0$ ). For every value of $d_0$ , an ensemble of eight independent simulations have been performed, starting from random initial perturbations to the flat film:

(4.1) \begin{equation} d(z,0) = d_0 + 10^{-3} \sum _{m=1}^{10} A_m\,\textrm{cos}(2 \pi m z R/\varLambda _{\textit{RP}} + \varphi _m), \end{equation}

where $A_m$ and $\varphi _m$ are uniformly distributed random numbers between $[0,1]$ and $[0, 2 \pi ]$ , respectively. Each simulation is run for a sufficiently long time so that the film attains its asymptotic profile. We use the one-way coupled WRIBL model (airflow does not affect the interface evolution, see § 3.2) and solve only the mucus-film equations (2.32)–(2.33), along with the constant cilia velocity boundary condition (2.12). Cilia transport, though, has no effect on the depleted zones because it does not alter the shape of the interface (see § 3.1). (We have checked that decreasing the cilia velocity by a factor of ten leaves the extent of depleted zones unchanged.)

In our simulations, the depleted zones are occupied by the precursor film. So, we measure the fraction of the wall occupied by depleted zones, $D_f$ , by measuring the length of wall along which $1-d\lt \delta$ , where $\delta \approx 0.01$ is the thickness of the precursor film.

Our results are presented for a range of initial film thicknesses $1-d_0$ in figure 8(a). Here, the markers and bars depict the mean and standard deviation of $D_f$ obtained from the ensemble of simulations for each value of $d_0$ . Surprisingly, adding more mucus (increasing $1-d_0$ ) does not decrease the extent of depleted zones, but rather increases it and leaves more of the wall exposed to particles.

Figure 8.

(a) Fraction of the wall occupied by mucus-depleted zones. The prediction (4.2) is compared with measurements from randomly initialised simulations. The markers and bars correspond to the mean and standard deviation of the values measured from the ensemble of eight simulations for each value of $d_0$ . (b–d) Shapes of unduloids (shaded in grey) for three values of the initial film thickness, overlaid with the final profile of the film (red-dashed line) from one of the runs in the corresponding ensemble of simulations. The same range of the axial $z$ coordinate is used in these three panels for ease of comparison; the cyan-shaded padding is added to $\varLambda _{\textit{RP}}/R$ to achieve the same $z$ -range.

This counter-intuitive variation of $D_f$ can be predicted by assuming that all the fluid contained in the initially flat film, of thickness $(1-d_0)R$ and length $\varLambda _{\textit{RP}}$ (the domain length), is ultimately gathered into an equilibrium-shaped unduloid. This unduloid’s axial width $\varLambda _U$ will be less than $\varLambda _{\textit{RP}}$ (both of which depend on $d_0$ ), and the difference $\varLambda _{\textit{RP}}-\varLambda _U$ will be the extent of the mucus-depleted zone. We thus obtain the following prediction for $D_f$ :

(4.2) \begin{equation} D_f \approx 1-\varLambda _U/\varLambda _{\textit{RP}}. \end{equation}

Here, $\varLambda _{\textit{RP}}$ , obtained from a linear stability analysis of the WRIBL model, is well approximated by the inviscid result of Rayleigh (Reference Rayleigh1892), $\varLambda _{\textit{RP}} = 2\pi 2^{1/2} d_0 R$ (Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015). To obtain $\varLambda _U$ , we calculate the shape of the unduloid $g(z)$ that (i) has a constant mean-curvature and (ii) encloses the entire volume of the film:

(4.3) \begin{align} \frac {1}{g}-\frac {(\textrm{d}_z g)^2}{2g}- \textrm{d}_{zz}g = K, \quad g(0) = 1, \quad \textrm{d}_z g |_0 = 0, \quad g(\varLambda _{U}/R) = 1, \nonumber \\ 2 \pi \int _{0}^{\varLambda _{U}/R} \int _{g(z)}^{1} r\, {\textrm{d}}r\, {\textrm{d}}z = \pi (1-d_0^2)\varLambda _{\textit{RP}}/R, \end{align}

where $K$ is a constant that is determined by the volume constraint. Note that we have used the $\mathcal{O}(\epsilon ^2)$ expression (2.7) for the mean curvature to be consistent with the WRIBL model. We solve (4.3) numerically using a shooting method. We start with a value for $K$ and solve an initial value problem for $g$ , stepping forward in $z$ until we satisfy the end boundary condition $g=1$ ; we thus obtain the unduloid length $\varLambda _U$ , after which the unduloid’s volume can be calculated. By varying $K$ , we obtain unduloids of all volumes less than $1.73 \pi$ , beyond which solutions do not exist. These unduloid shapes match quite well with those computed by Everett & Haynes (Reference Everett and Haynes1972) using the full expression for the mean curvature; importantly, the dependence of $\varLambda _U$ on the initial thickness of the film ( $1-d_0$ ) is well approximated (see the supplementary material). This is not the case if one uses the linearised expression for the curvature (obtained in the limit $1-d_0 \ll 1$ ), as is done in traditional lubrication models of thin films: the resulting equation is linear in $g$ and yields unduloids with a fixed, $d_0$ -independent, width of $2\pi$ (Lister et al. Reference Lister, Rallison, King, Cummings and Jensen2006). So, in addition to pinchoff dynamics (Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015), the shape of unduloids is another qualitative feature that necessitates the use of the $\mathcal{O}(\epsilon ^2)$ WRIBL model.

The variation of $D_f$ predicted by (4.2) is plotted as a solid-black line in figure 8(a) and seen to agree rather well with the simulations. Importantly, the increase in $D_f$ with the mucus volume can now be understood in terms of the change in the unduloid’s shape, which is illustrated in figure 8(b–d). Here, we also overlay the final profile of the film (red-dashed line) from the simulations, on top of the unduloid (grey shading); despite the use of a precursor film, the simulated interface profile agrees rather well with the conceptual ideal of an equilibrium unduloid and a fully depleted zone. Now, as the film thickness increases, destabilising capillary forces due to the radial curvature increase in strength relative to stabilising axial-curvature forces. As a result, the unduloids (and the corresponding humps in the simulations) become deeper, so much so that their length $\varLambda _U$ reduces even as the total volume increases. The wavelength of the fastest-mode $\varLambda _{\textit{RP}}$ also decreases, but not as much as $\varLambda _U$ . The net effect therefore is that a larger portion of the wall is left devoid of mucus.

The small variation in $D_f$ among the different random simulations (see the small error bars in figure 8 a) confirms our expectation that regardless of the details of the initial perturbation (so long as it triggers the fastest-growing mode of the Rayleigh–Plateau instability), the film ultimately evolves into an unduloid-shaped hump with an adjacent depleted zone. A complete understanding of the role of random perturbations, though, requires an analysis of simulations on long domains that span multiple wavelengths of the fastest growing mode. Our ongoing work in this direction shows that (4.2) remains a good approximation to $D_f$ even on long domains. We comment further on this point in our concluding remarks in § 6.

An important consequence of (4.2) is that the minimum value of $D_f$ , attained in the limit of a vanishingly thin film ( $d_0 \rightarrow 1$ ), is far from zero – approximately 0.30 as seen from the solid line in figure 8(a). (This limiting value is what is predicted when using the linearised expression for the curvature, as shown in the supplementary material.) So, an appreciable fraction of the wall is always depleted of mucus, leaving it vulnerable to airborne particles. Additionally, adding more mucus only worsens the situation by producing deeper but shorter humps that expose even more of the wall. Does this increased exposure imply that more particles will deposit on the wall when more mucus is present, or will the deeper humps be able to trap more particles? We directly address this question, for particles spanning a wide range of sizes, in the next section.

5. Particle transport and deposition

We consider two different scenarios in this study of particle deposition, which correspond to two extremes of inter-breath particle turnover. In the first scenario, the same particles persist within the airway throughout the simulation, which spans 30 breathing cycles (30 s). In the second scenario, we assume that all particles are replaced by fresh ones after each breath, due to mixing between tidal air and reserve air or between neighbouring airways (Altshuler et al. Reference Altshuler, Palmes, Yarmus and Nelson1959; Wang Reference Wang2005, chap. 5).

A number concentration of $300\,\unicode{x03BC}\textrm{g m}^{-3}$ , which is typical for $\textrm{PM}_{2.5}$ (epa.gov 2024), translates to approximately 16 micron-sized particles in the airway under consideration. So, in the first scenario, we introduce 16 particles into the periodic airway domain, at random locations, and track their motion for 30 s (while removing particles that deposit on the film). We then perform an ensemble of 125 such simulations, each starting from a different random initial distribution of particles. We thus obtain good statistics for deposition results, such as $\phi (t)$ , the fraction of particles that have deposited (on the wall or mucus) at any time during the simulation. Now, because the particles do not interact with each other in the simulation, the mean value of $\phi$ obtained from the ensemble of $125$ simulations would be identical to the value obtained from a single simulation with $16\times 125 = 2000$ particles. Hence, the chosen number concentration, i.e. the number of particles introduced into the airway in each simulation, will only affect the variance in the deposition fraction and not the ensemble-averaged mean results.

In the second scenario, we again introduce 16 particles at random locations in the airway and track them for one breathing cycle; we then remove all non-trapped particles and introduce a fresh set of 16 particles. This process is repeated for 30 breathing cycles (30 s of respiration). We then perform an ensemble of 60 such simulations. Of course, because the particles are replaced after each cycle, this calculation is equivalent to running an ensemble of 1800 simulations, each for one breathing cycle.

From the perspective of particle transport, the key difference between the two scenarios is that only in the first scenario will the distribution of particles have sufficient time to respond to the non-homogeneous nature of the flow field. For example, the curved streamlines of the airflow near the mucus hump will induce a slight centrifugal drift of weakly inertial particles towards the centreline, which in the first scenario will add up from breath to breath to produce an accumulation of particles near the centreline. Such accumulative effects will not occur in the second scenario, wherein complete turnover of particles takes place after each breath.

While discussing our results, we shall distinguish between deposition on the mucus humps (mucus entrapment) and deposition on the mucus-depleted zones (wall deposition). In the latter case, the particles deposit not directly on the epithelial cells of the airway wall, but rather on the PCL layer which our model ignores. Nevertheless, we shall, for simplicity, use the term wall deposition because these particles are much more likely to reach the wall through the thin, low-viscosity, PCL layer than through the combination of the thick, viscous, mucus hump followed by the PCL layer.

In § 3.2, we had noted that the air flow was orders of magnitude faster than the cilia-induced translation of the interface and so we had anticipated that ciliary transport will not affect particle deposition. We have checked this and indeed found that our deposition results remain unaffected by decreasing $\langle u_c \rangle$ by a factor of ten. So, though ciliary transport is present in the particle-transport simulations, it does not contribute to deposition.

Let us begin with the first scenario and examine how the deposition fraction $\phi$ depends on the particle size. Results obtained after 30 breathing cycles (30 s) are presented in figure 9(a) for a relatively thick film. Three curves of $\phi$ are shown, corresponding to the fraction of particles deposited on the mucus hump, on the airway wall (depleted zone) and on either surface (total deposition). The particle diameter $d_p$ is represented by $\textit{Pe}$ and $St$ (see (2.38) and table 3); recall that $\textit{Pe}$ represents the importance of convection relative to Brownian diffusion, while $St$ indicates the importance of particle inertia. Both $\textit{Pe}$ and $St$ increase with size ( $Pe \sim d_p$ and $St \sim d_p^2$ ), so that the particle motion changes from diffusive to tracer-like to inertial. This varying nature of particle motion may be appreciated by comparing supplementary movie 1 with movie 2 and movie 3, which illustrate the motion of particles with $Pe = 2 \times 10^4$ ( $St = 9 \times 10^{-7}$ ), $Pe = 9 \times 10^5$ ( $St = 2 \times 10^{-3}$ ) and $Pe = 9 \times 10^6$ ( $St = 2 \times 10^{-1}$ ), respectively. (These movies use a large number of particles to better visualise the qualitative features of their motion.)

Figure 9.

Deposition fraction $\phi$ as a function of particle size, after 30 breathing cycles (scenario one), in an airway with (a) a relatively thick film, $1-d_0=0.115$ , and (b) a very thin film, $1-d_0 = 0.02$ . The particle size is represented by the corresponding values of $St$ (bottom axis) and $\textit{Pe}$ (top axis). Snapshots of the two films, in their long-term asymptotic states, are shown as insets in the corresponding panels. For each particle size, the marker represents the mean of the distribution of $\phi$ , obtained from the ensemble of 125 simulations; the bar corresponds to the interquartile range which is indicative of the spread of the distribution. The fraction of particles depositing on the mucus hump, the depleted zone of the wall and on either surface are shown separately (see the legend).

Figure 9(a) shows that the total deposition initially decreases with particle size, as the strength of Brownian forces reduce. However, there is an increase of deposition for the larger, inertial particles. That this significant increase is associated with the presence of a deep mucus hump is demonstrated by comparing figures 9(a) and 9(b). In the latter, which corresponds to a much thinner film with a shallow hump, the deposition fraction remains very small even for inertial particles.

The non-monotonic variation of total deposition with particle size, seen in figure 9(a), is a well-known feature of experimental measurements of deposition in the lungs (Morawska et al. Reference Morawska, Buonanno, Mikszewski and Stabile2022). However, the increase for inertial particles is much stronger in entire-lung measurements because the effects of particle inertia are very strong in the turbulent flows of the upper airways, as well as in the secondary circulatory flows at airway bifurcations (Guha Reference Guha2008). Here, we see that a non-monotonic variation of deposition, albeit weaker, arises even in a single airway provided the mucus film is thick enough to form prominent humps.

Figure 10.

Evolution of the ensemble-averaged deposition fraction on (a) the depleted zones of the wall and (b) the mucus hump. Results are shown for six different particles sizes; for convenience, both the $\textit{Pe}$ and $St$ values of these particles are mentioned in the legend, which is split across the two panels. The solid-black line corresponds to a growth of $t^{1/2}$ .

The variation of the deposition fraction with time, for different particle sizes, is presented in figure 10. The small particles, whose deposition is driven by Brownian motion, exhibit a $t^{1/2}$ growth of $\phi$ at early times (this scaling can be derived by considering diffusive mass loss into an absorbing circular boundary as outlined in the supplementary material). Given sufficient time, all these particles will be deposited on either the wall or the mucus humps. The effects of inertia are first seen for $St = 2 \times 10^{-3}$ , for which $\phi$ saturates to a low value (figure 10 b). As shown later, the cessation of further deposition is a result of gradual centrifugal migration towards the centreline, owing to the curved air streamlines below the mucus hump. As these inertial forces become stronger, the particles veer-off from streamlines more violently and collide with the mucus hump; hence, the increase in deposition for the largest particle with $St = 2 \times 10^{-1}$ (figure 10 b).

Returning to figure 9(a), let us examine the difference between the deposition on the mucus hump (blue dashed line) and on the exposed wall (red solid line). For the small particles, we see that though a higher fraction deposits on the mucus hump, a significant fraction does reach the wall. The difference between the deposition fractions mirrors the difference in the surface areas of the mucus hump and the exposed wall (55 % and 45 %, respectively, of the total surface area exposed to air). This correspondence is a sign of the diffusive nature of deposition of the small particles (see in particular the values for the smallest particle with $Pe = 2 \times 10^2$ ). The deposition of inertial particles on the hump and the wall is governed by other physical effects and will be discussed later in this section.

Figure 11.

Effect of the initial film thickness (mucus volume) on the deposition fraction. The results from figure 9(a,b) are overlapped for deposition on (a) the depleted zone of the wall and (b) the mucus hump.

Next, let us examine the effect of the mucus volume fraction with the aid of figure 11, which overlaps the results for the two film thickness shown previously in figure 9(a,b). We have already seen that thicker films produce larger mucus-depleted zones (see figure 8 a). One would expect the widening of depleted zones to allow more small Brownian particles to deposit on the wall when the film is thicker. This is exactly what we observe in figure 11(a). Consequently, a smaller fraction of small particles is entrapped by the mucus humps in the case of the thicker film (figure 11 b). These differences between the two films are less for the larger tracer-like particles, which diffuse more slowly. However, deposition of the largest, inertial particles is again markedly higher for the thicker film – especially on the mucus hump – owing to inertial effects.

Inertia causes the large particles to deviate from the air streamlines and collide with the wall and the mucus hump. When heavy particles approach the mucus hump, their inertia prevents them from going around the hump along with the curved streamlines. Instead, they continue moving straight and collide with the mucus hump. Evidence for such inertial impaction (Guha Reference Guha2008) is provided by the collision velocity, $V_{||}$ , which is the relative velocity between the particle and the deposition site (mucus hump or wall) at the time of deposition, projected along the approach direction. Figure 12 presents the mean collision velocity of particles with the mucus hump (panel a) and the wall (panel b) for the case of the thick film (figure 9 a); here, the bars represent the inter-quartile range which is indicative of the spread of the distribution of $V_{||}$ . The smallest diffusive particles have a large $V_{||}$ because of the dominance of fluctuating Brownian forces. The tracer-like particles have a very low $V_{||}$ , as expected for particles whose velocity is almost entirely aligned with the air streamlines which are tangential to the mucus hump and the wall. The relatively large, inertial particles though show a marked increase in $V_{||}$ , especially for collisions with the mucus hump which protrudes into the airway (figure 12 a); as anticipated, inertial particles deviate from the streamlines and collide with the mucus hump.

The slight decrease in $V_{||}$ for the collisions of the largest particle with the hump (compare the second-largest with the largest particle in figure 12 a) signifies the growing influence of the particle interception effect (Guha Reference Guha2008). Because of its increasing radius, the centre of the particle need not deviate as strongly from the airflow streamlines for the particle surface to make contact with the mucus film. Thus, deposition will include more low-velocity collisions as the particle size increases.

Figure 12.

Collision velocity $V_{||}$ of particles with the depositing surface: (a) mucus hump; (b) depleted-zone of the wall. The marker and the bar represents the mean and the inter-quartile range (which measures the variation about the mean) of the distribution of $V_{||}$ for each particle size. These results correspond to the case of the relatively thick film in figure 9(a).

The collision velocity of inertial particles with the wall is much smaller than that with the hump, even for the largest particle (compare figure 12 a,b). Deposition on the wall, though, is aided by the non-uniform number concentration of inertial particles in the airway: a subtle interplay of varying airflow speed and streamline curvature leads to an elevated number concentration near the depleted zone on the wall.

The spatial distribution of particles is examined in figure 13. The mean radial location of all non-trapped particles, $\langle r \rangle _p$ , is plotted as a function of time in figure 13(a). A gradual migration towards the centreline, $r = 0$ , is apparent here for the relatively large inertial particles; such migration does not occur for the small diffusive particles. The small oscillations in $\langle r \rangle _p$ , present for all particles, are a result of the axisymmetric geometry and the converging–diverging nature of the air streamlines, which causes a small convergence towards the centreline, followed by a corresponding divergence, each time particles are carried past the hump. These oscillations do not lead to a net convergence, though, and the preferential concentration of inertial particles near the centreline is caused by gradual centrifugation. Particles that start out at intermediate radial locations – neither near the centreline nor near enough to the mucus film to be intercepted by the mucus hump – will be transported to and fro across the airway by the respiratory air flow. As these particles pass below the mucus hump, the curved nature of the streamlines (see the inset of figure 9 a) will result in the particles being centrifuged towards the centreline. This migration occurs slowly, though, and so it takes several breathing cycles for the particles to accumulate at the centreline (figure 13 a).

Figure 13.

Non-homogeneous distribution of particles in the airway, for the case of the relatively thick film (see figure 9 a). (a) Evolution of the mean radial location of particles (non-deposited) in the airway; the migration of inertial particles ( $St = 0.009,\,0.2$ ) towards the centreline is clearly visible. (b) Probability distribution function of the radial location of particles, integrated over time (data from $t = 0$ to $t = 30T_b$ is used to construct the distribution). The distribution near the centreline is shown in the main panel, while that near the wall is show in the inset.

The probability distribution function (p.d.f.) of the radial location of non-trapped particles, integrated over time, is show in figure 13(b). The main panel reinforces the message of figure 13(a), that small particles are uniformly distributed, while inertial particles accumulate near the centreline. The inset, though, reveals a surprising accumulation of inertial particles near the depleted zone of the wall (there is no accumulation below the hump which extends up to $r \approx 0.6$ ). Supplementary movie 3 shows how heavy particles that begin near the depleted zone tend to remain there (though the majority of particles, which at some time pass by the hump, are centrifuged towards the centreline). This is due, in part, to the low air velocity near the depleted zone (which has the highest cross-sectional flow area). Over a breathing cycle, the particles travel only a short way past the edge of the hump before returning to the depletion zone, and hence they do not pass the bottom of the hump where the centrifugation towards the centreline acts most strongly. In contrast, because the curvature of the streamlines near the edge of the depletion zone (where it meets the hump) is opposite to that of the streamlines below the hump, the inertial particles are weakly pushed towards the edge of the depletion zone. The extent of this accumulation near the wall is much weaker than that near the centreline (evidenced by figure 13 b) because the streamline curvature and the air flow velocity is weaker near the wall than below the mucus hump.

To summarise the results of the first scenario, there is a non-monotonic variation of deposition with particle size because small particles deposit due to Brownian forces and large particles deposit due to inertial effects. The latter are rather subtle: strong centrifugal forces near mucus humps produce high-velocity collisions, whereas weak centrifugal forces near mucus-depleted zones of the wall cause particles to accumulate and then gradually deposit there.

Figure 14.

Variation of the deposition fraction $\phi$ with particle size, for the second scenario in which particles are replaced after each breath. Results for the thick and thin film are compared for deposition on (a) the mucus hump and (b) the wall. The insets are zooms that show the variation of $\phi$ for the large particles. The markers and the bars represent the mean and the inter-quartile range of the deposition fraction measured from an ensemble of 60 simulations, each run for 30 breathing cycles (a total of 1800 breaths between which particles are entirely replaced).

Let us now turn to the second scenario, in which undeposited particles are removed at the end of each breathing cycle and a new set of 16 particles are introduced at random locations. First, it is important to note that the manner in which $\phi$ is calculated in the second scenario will typically result in smaller values as compared with the first scenario. For the sake of illustration, let us assume that the same number of particles deposit onto the hump in each cycle, for both scenarios. Even so, the second scenario will yield a much smaller deposition fraction because the fraction $\phi$ is calculated by dividing the total number of deposited particles by $30\times 16$ and not by $16$ (as would be done for the first scenario one). In this context, a fair way to compare deposition between the two scenarios is to examine the actual number of particles deposited onto the hump after 30 breathing cycles, which is given by $16 \phi$ and $30*16 \phi$ for scenarios one and two, respectively. Note that if the deposition is rapid, such that nearly all the particles introduced into the airway deposit in a single cycle, then the two scenarios will yield comparable values of $\phi$ ; we shall find this to be the case for the smallest diffusive particles. In the following discussion of the second scenario, we shall examine how $\phi$ changes with the particle size, in light of the corresponding variation in the first scenario.

Because the particles are replaced very frequently, the inertial particles will not have time to accumulate near the centreline or near the wall, as they did in scenario one (figure 13 b). Therefore, while particle inertia should still promote deposition on the mucus hump due to centrifugation and high-velocity collisions, there should not be any wall deposition due to inertia because particles will not accumulate near the wall in scenario two. This is exactly what we see in figure 14, which compares the results for deposition on the mucus hump (panel a) and on the mucus-depleted portion of the wall (panel b), for the two cases of a thick and a thin film (see the legend). More inertial particles are trapped by the thicker film with the deeper mucus hump (see the inset of figure 14 a and note the significant relative increase of $\phi$ for the largest particle). As anticipated, there is no such increase in the case of wall deposition of inertial particles (see the inset of figure 14 b). For small diffusive particles, we see that, just as in scenario one, the thicker mucus film allows more particles to deposit on the wall because of the presence of larger mucus-depleted zones (figure 14 b).

The key takeaway from the results of the second scenario is that particle turnover counter-acts the effects of inertia-driven preferential accumulation. So, if a substantial fraction of particles in an airway are replaced after each breath (thanks to mixing between airways or between tidal and reserve air volumes), then a higher mucus volume will not cause more wall deposition of heavy inertial particles. Taken together, the two scenarios suggest that the primary effects of increasing the mucus volume fraction are more mucus entrapment of large inertial particles and more wall deposition of small diffusive particles. Which of these opposing outcomes is detrimental or beneficial to human health depends on the nature of the particles – harmful allergens and pathogens, or therapeutic drugs.