## 3 Results and Discussion
### Parametric Investigations
Given that the reported value for the spacing between cilia is 0.35 \(\upmu\)m [59], further decreasing the spacing between two adjacent cilia can cause spurious numerical effects since the neighbouring cilia may overlap, and as a result, two Lagrangian points with opposite velocities can occupy the same fluid grid cell. It is also important to emphasise that the wavelength formed by motile cilia is kept constant for all of the ciliary densities. Thus, the phase numbers of cilia are interpolated according to their positions within the ciliary domain. The linear dependency of \(\langle Q\rangle\) on frequency (Figure 8a) and number of cilia (Figure 8e) can be expected as the driving force linearly increases with frequency and number of cilia in the creep flow regime and therefore the average PCL flow would linearly depend on  
Figure 9: (**a**) Schematic of ciliary pack. (Pack 1): equilateral triangular; (Pack 2): isosceles triangle; (Pack 3): square. (**b**) Schematic of an array of 144 cilia with three patches representing the regions of missing cilia. these two parameters. Figure 8f shows the relation between the ciliary density and \(\langle Q\rangle\). This simulation results show that although there is a decrease for the case with 64 cilia, the flow rate initially increases with increasing the number of cilia because the total driving force increases. However, when the cilia are too close to each other the flow generated by one cilium would be partially cancelled out by the surrounding cilia and one can therefore see a plateau in \(\langle Q\rangle\) as the number of cilia approaches 400. In all of the simulations presented, results obtained for the 2 nd and 3 rd cycles have been almost identical, regardless of packing type, patch size, and ciliary density. Hence, the time-averaged results presented in this section are obtained over the cycle 3 with the period of \(T=0.1\) s.  
It is now interesting to investigate the variation of time-averaged surface mean velocity \(\langle u_{surf}\rangle\) for all of the six above-mentioned parameters (see Figure 10). One can clearly see that the curves exhibit relatively similar trends as the time-averaged flow rates \(\langle Q\rangle\), particularly in parameters (I), (IV), (V), and (VI); Figure 10a,d-f. This similarity is to be expected since the flow rate is the product of area and velocity. The linear relationship between \(\langle u_{surf}\rangle\) and \(CBF\) in Figure 10a is also to be expected due to the negligible inertia forces (very low Reynolds number) for this problem which, consequently, results in a linear momentum Equation (1). Additionally, the experimental findings also indicate that ciliary forces depend linearly on the \(CBF\)[60]. Figure 10b shows that the surface mean velocity depends highly on the ciliary length \(L\). Decreasing \(L\) from 7 to 4 \(\upmu\)m, for example, could reduce \(\langle u_{surf}\rangle\) by 206%, indicating that smoke-induced shortening of cilia could decrease the muco-ciliary clearance significantly. Interestingly, it was found that the standard beat is able to slightly improve the surface mean velocity compared with the planar beat (see Figure 10c), although it had less fluid transport capability. On the other hand, \(\langle u_{surf}\rangle\) experienced a reduction of over 80% when using the stiff planar beat, and the surface remained almost stationary for the wiper motion. As mentioned in above, \(\langle u_{surf}\rangle\) exhibits relatively similar trend as \(\langle Q\rangle\) when comparing different ciliary pack(parameter (IV)). That is to say, Pack 1 has the largest value, followed by Pack 2 and Pack 3, respectively (see Figure 10d). [PERSON] et al. [61] measured the in-plane fluid velocities above the epithelium surface of mucus-free mice tracheae using the micro Particle Image Velocimetry technique and observed that regions where no beating cilia are located had smaller velocities compared to the regions with ciliated cells.