## 2 Numerical Method
The staggered marker-and-cell (MAC) Eulerian mesh [42] is used for the spatial discretization (as shown in Figure 2c), and the intermediate velocity field \(\vec{\mathbf{u}}^{*}\) is calculated as follows:  
\[\
abla^{2}\vec{\mathbf{u}}^{*}-\lambda\vec{\mathbf{u}}^{*}=RHS, \tag{7}\]  
where \(\lambda=2\rho/\mu.\Delta t\) and  
\[RHS=\lambda\vec{\mathbf{u}}^{n}-\
abla^{2}\vec{\mathbf{u}}^{n}+\frac{2}{\mu} \big{(}\rho(\vec{\mathbf{u}}\cdot\
abla\vec{\mathbf{u}})^{n}+\
abla p-\overrightarrow {F}^{n+1}\big{)}. \tag{8}\]  
Since the \(\vec{\mathbf{u}}^{*}(\sim\vec{\mathbf{u}}+\
abla\phi)\) does not satisfy the divergence free condition in Equation (2), a pressure increment \(\phi^{n+1}\) is computed as  
\[\
abla^{2}\phi^{n+1}=\frac{\
abla\vec{\mathbf{u}}^{*}}{\Delta t},\ \ \mathbf{n} \cdot\
abla\phi^{n+1}|_{\partial\Omega}=0. \tag{9}\]  
Figure 2: (**a**) Flow domain with an immersed boundary (IB). \(\Delta S\) is the spacing between the control points, and \(h=\Delta x=\Delta y=\Delta z\) is the grid spacing. (**b**) Interpolation strategy of the immersed boundary method. (**c**) The MAC mesh cell. The velocity components \(u\), \(v\), and \(w\) are defined at the cell faces, while the pressure \(p\) is introduced at the cubic cell centre. To calculate the Helmholtz equation for \(\vec{\mathbf{u}}^{*}\) in (7) and the Poisson equation for \(\phi^{n+1}\) in (9), we employ FISHPACK [43] fast solvers. Having solved for the \(\vec{\mathbf{u}}^{*}\) and \(\phi^{n+1}\), the velocity and pressure fields are then updated as  
\[\vec{\mathbf{u}}^{n+1}=\vec{\mathbf{u}}^{*}-\Delta t\
abla\phi^{n+1}, \tag{10}\]  
\[p^{n+1/2}=p^{n-1/2}+\phi^{n+1}-\frac{\mu}{2\rho}\
abla\vec{\mathbf{u}}^{*}. \tag{11}\]  
The standard projection method has been successfully applied in many muco-ciliary studies [6; 8; 23; 44; 45], and all the authors have reported satisfactory values for the fluid transport when compared with the available experimental data. Nonetheless, the accuracy of this method at low Reynolds numbers has been questioned because the projection on divergence-free fields has been observed to induce tangential spurious velocities [46]. Although this is an important limitation to consider, the actual degree of loss of accuracy that might result is unknown for muco-ciliary systems. What makes the projection method attractive is that, at each time step, one only needs to solve a sequence of decoupled elliptic equations for the velocity and the pressure, which makes it more computationally efficient. However, it is necessary to evaluate the accuracy of the algorithm which is tested for the Couette flow, as well as for Stokes flow, in a lid-driven cavity, for \(Re=0.01\), which is a creeping flow condition of analogous flow regime to muco-ciliary transport. Figure 3a provides results for the Couette flow at a range of times compared with the analytical solution of the Navier-Stokes equations, which can be expressed as [47]  
\[u(y,t)=\frac{u_{0}}{L_{y}}y+\sum_{n=0}^{\infty}\frac{2u_{0}}{n\pi}(-1)^{n}sin \bigg{(}\frac{n\pi}{L_{y}}y\bigg{)}\exp\bigg{(}-\
u\frac{n^{2}\pi^{2}}{{L_{y} }^{2}}t\bigg{)}, \tag{12}\]  
where \(L_{y}\) is the width of the channel, \(\
u\) is the fluid kinematic viscosity, and \(u(y,t)\) is the fluid velocity calculated at position \(y\) and time \(t\). Comparable to the results obtained by the analytical solution, the present results deviated by no more than \(0.12\%\) which is acceptable within the current context. As for the second validation (Figure 3b), the streamline results of the present work are in good agreement with those obtained by an analytical method in the works of [PERSON] [48] and [PERSON] [49], supporting the use of the numerical algorithm to simulate Stokes flows for our purposes. Detailed discussions on the validity and robustness of the implemented IB method are found in Reference [23; 44; 50] and hence not repeated here. Table 2 summarizes the computational steps used in the present work.