## 3 Results and Discussion  
### Cilia-Driven Fluid Flow  
To assess the numerical convergence of the discretization scheme, we evaluate the effects of time step, fluid grid, and number of control points on the time-averaged fluid velocity in the \(x\)-direction. This is computed over the total PLC volume using \(\langle u\rangle=1/(TV_{f})\iint u(x,y,z,t)dV_{f}dt\) where \(V_{f}\) is the fluid volume and \(T=0.3\) s. The results show that the solution changes little under spatial and temporal mesh refinement and the discretization scheme reaches convergence over the third cycle, as seen in Figure 4. Therefore, the fluid grid size \(h=0.5\) (where \(h=L_{\Omega}/m_{x}=W_{\Omega}/m_{y}=H_{\Omega}/m_{z}\)), the control points \(N=20\) and the time step \(\Delta t=1\times 10^{-6}\) s are chosen for our simulations. Figure 5 shows the time-averaged velocity \(\langle\mathbf{u}\rangle=(1/T)\int\mathbf{u}dt\) generated by ciliary motion. The mean velocity contour in the middle section of the computational domain in the \(z\)-direction reveals a very complex region, as depicted in Figure 5a. In addition to this, since the cilia beat most effectively in the \(x\)-direction, it is expected to have regions with higher velocities in this direction. The largest values of the mean velocity magnitude \(\langle\mathbf{u}\rangle\) are observed on the \(z=8\) um plane (Figure 5b) which can facilitate the mobility of the ML. Figure 5c,d depict the \(x\)- and \(y\)-component of velocity field (i.e., \(\langle u\rangle\) and \(\langle v\rangle\)) at two representative planes: \(y=0\) and \(x=0\), respectively. It is observed that the maximum contour value of \(\langle u\rangle\) is approximately four times larger than that of \(\langle v\rangle\). To quantify these flow changes, the time-averaged velocity profiles as functions of \(z\) are constructed in the middle of the computational domain (intersection of \((x,z)\) and \((y,z)\) planes), as seen in Figure 5e. These findings clearly confirm that a considerable amount of flow in the direction parallel to effective stroke is generated. Referring to Figure 5c, one can also observe that the main direction of the fluid transport above the cilia tips is towards the positive \(x\)-direction which agrees well with the actual clearance mechanism; unidirectional transport of the ASL over time. Furthermore, our calculations for the time-averaged fluid velocity which is computed over the total PLC volume in \(y\)-direction \(\langle v\rangle=1/(TV_{f})\iint v(x,y,z,t)dV_{f}dt\) led to  
Figure 4: Numerical convergence study by varying the time step size \(\Delta t\), number of fluid grids (i.e., \(m_{x},m_{y},m_{z}\)), and number of control points \(N\). a very small value: \(\langle v\rangle=0.145\)\(\upmu\)m/s vs. \(\langle u\rangle=3.68\)\(\upmu\)m/s. Henceforth, we shall present our findings for flow parallel to the direction of cilia effective stroke. Figure 6 shows streamlines and velocity profiles generated by ciliary actions during the stable beating cycle which takes 0.1 s. The streamlines are taken at the plane \(z=8\)\(\upmu\)m (left panel) and at the plane \(y=1.5\)\(\upmu\)m (shown by a dotted line in the top view). Six equally spaced velocity profiles are  
Figure 5: Time-averaged flow velocity. (**a,b**) Contour plots of mean velocity fields at planes \(z=4\)\(\upmu\)m and \(z=8\)\(\upmu\)m, respectively. (**c,d**) Velocity contours at planes in the directions parallel (\(y=0\)\(\upmu\)m) and perpendicular (\(x=0\)\(\upmu\)m) to effective stroke, respectively. (**e**) Flow profiles in the \(x\)- and \(y\)-directions as functions of \(z\). The profiles are plotted at the intersection of planes shown in (**c,d**) of this figure. shown with vertical dotted line for \(-24\leq x\leq 24\) um (right panel). The cilia are labelled 1-6 from left to right. We now seek to understand the interaction between the fluid and cilia at various time instants.