## 3 Results and Discussion
### Cilia-Driven Fluid Flow
Viewing from above, many small size vortices are formed above the cilia tips due to the packed distribution of cilia, as well as the complex nature of cilia-driven fluid field. Such turbulence-like streamlines are due to the flow separations in the vicinity of ciliary actions. However, the overall fluid regime is still laminar. Outside the ciliary domain (i.e., \(|(x,y)|\geq 16.5\) um), it is more likely to have bigger recirculation zones. The elevation view shows a row of six cilia which are at different beating phase, depending on the simulation time. Case (a) shows a snapshot of the results at the beginning of the beating cycle, i.e., \(t=0\) s. In case (b), cilium 1 is seen to stand up and perform its effective stroke. The action of this cilium is pushing the fluid forward in the positive \(x\)-direction; the corresponding velocity profile shows this positive fluid motion. The cilium continues to move forward (case (c)), but compared to the previous case, the corresponding fluid velocity profile shows smaller values. Case (d) depicts the recovery stroke for this cilium where an opposing vortex is set up in the fluid above the cilium tip, resulting in a reverse flow with a negative velocity profile. One can note that this relatively strong vortex is formed above the other cilia with similar beating phase in the other snapshots. Eventually, at \(t=0.08\) s, which corresponds to the case (e), the cilium 1 start reaching to its previous initial position to begin the next cycle. This biological fluid-structure interaction is very similar to what is observed in nature. For example, birds extend their wings to overcome the viscous resistance and then bend them closer to their body to avoid cancelling out their effective forward motion. To provide information on the PCL transport, the flow rate is computed as  
\[Q=\iint u(x=0\mu m,y,z,t)dydz. \tag{13}\]  
The surface mean velocity where the ML sits on top is also obtained as  
\[u_{surf}=(1/A)\iint u(x,y,z=8\mu m,t)dA. \tag{14}\]  
The fluid-flow induced wall shear stress at the epithelia are also another important factor of interest that, to our knowledge, has not been numerically investigated, which is calculated as  
\[\tau_{w}=(\mu/A)\iint\ \left((\partial u(x,y,z=0)/\partial z)^{2}+(\partial v(x,y,z=0)/\partial z)^{2}\right)^{1/2}dA. \tag{15}\]  
It has been experimentally observed that the shear stress can alter the epithelial paracellular permeability (i.e., diffusion of ions or molecules via the intercellular space) significantly [53]. The magnitude of fluid vorticity is the last metric of interest which is averaged over the PCL volume and computed as  
\[|\omega|=(1/V_{f})\bigg{|}\iiint(\
abla\times\mathbf{u})dV_{f}\bigg{|}. \tag{16}\]  
Figure 7 shows the transients of flow rate, PCL surface average velocity, wall shear stress, and volume-averaged vorticity magnitude obtained from Equations (13)-(16). With 72 cilia beating at 10 Hz, the average of PCL flow rate, surface mean velocity, wall shear stress, and vorticity magnitude are obtained as follows: \(Q=1533\) um\({}^{3}\)/s, \(u_{surf}=8.19\) um/s, \(\tau_{w}=0.08\) dynes/cm\({}^{2}\), and \(|\omega|=0.95\) s\({}^{-1}\), respectively. Our surface velocity results are qualitatively consistent with those reported in Reference [54], but different quantitatively, with the time-averaged PCL mean velocity being 3.18 um/s. The reasons for this discrepancy are due to differences in the adopted ciliary density and prescribed beating pattern. In addition, they considered a synchronized cilia beating (no methanol wave) and a fluid with variable viscosity, both of which can contribute to the further difference. Figure 6: _Cont._