### Longitudinal motion of Corti fluid

A 6 mm-long and 25 µm-thick fluid domain was subjected to traveling waves along one longitudinal boundary (Fig. 2A). For the case shown in panels (A–C), the wavelength was 4 mm and wave speed (phase velocity) was 80 m/s to the right (the longitudinal direction). The vibration amplitude was 25 nm. The color contour in panel A shows the pressure field when the peak inward (transverse) displacement occurred at the center of the tube. The streamlines (solid curves) depart from and arrive at the stimulating boundary. Velocity profiles at the midsection (vertical broken line in panel A) reveal that the flow moves back and forth in both directions as pressure waves propagate (panels B and C). The velocity in the longitudinal direction is greater than the transverse direction by two orders of magnitude for the simulated case. This implies that the fluid motion is primarily longitudinal. The longitudinal velocity that is set to zero at the stimulation boundary grows steeply over 5 µm to reach 300 mm/s. That is, approximately 75% of the fluid volume is under the influence of strong longitudinal motion. In contrast, the transverse velocity is greatest along the stimulating boundary and decreases toward the line of symmetry.

Long waves (*λ*/*r*_{0} > > 1) cause large longitudinal motion, as shown in Fig. 2D. In contrast, as the wavelength gets shorter (*λ*/*r*_{0} < 1), longitudinal and transverse velocities become comparable. For a long wave, there was minimal pressure gradient in the transverse direction, thus the fluid was forced to move longitudinally. On the contrary, for a short wave (*λ*/*r*_{0} < 1), the pressure varied appreciably in both transverse and longitudinal directions (results not shown). Mammalian cochleae have *r*_{0} < 50 µm and experience wavelengths ranging from 100 µm to a few millimeters, so the long-wave case shown in panels A–C represents the natural conditions better than the short-wave cases.

Overall fluid velocity was proportional to the vibration amplitude *a* and the wave speed *c* (= *fλ*). As shown in panel E, as long as the wave velocity is the same, combinations of different frequencies and wavelengths resulted in comparable fluid velocity. In the plot, fluid velocity was represented by the root-mean-square velocity over the entire fluid domain.

### Fluid trajectories obtained from Stokes drift velocity

The result in Fig. 2 shows large instantaneous longitudinal fluid velocity when long waves propagate at high speed. However, this does not necessarily result in longitudinal mass transport because fluid particle motions are primarily cyclic—they return to nearly the same position after each cycle. Effective long-term mass transport is caused by the offsets particles experience over many cycles, due to Stokes drift. For the case of Fig. 2A–C, the velocity field was analyzed to quantify those offsets (Fig. 3).

Figure 3A shows the paths of nine particles over one cycle (the sub-domain of panel A corresponds to the broken rectangle in panel C). As expected, particles return to nearly the same position after a cycle but experience small offsets. Corresponding Stokes drift velocity was between 0.1 and 1 mm/s over half of the fluid domain in the case of c = 80 m/s (Fig. 3C). This Stokes drift velocity is smaller than the peak instantaneous velocity by more than two orders of magnitude (cf. ~ 300 mm/s in Fig. 2B). Integration of the Stokes drift velocity **u**_{D} results in the trajectories in Fig. 3C–E (red curves) that approximate fluid particle motions after excluding fast oscillatory motion. The oval shape of particle paths is aligned with their major axis parallel to the longitudinal direction, consistent with the velocity ratio between the two directions. Stokes drift is a nonlinear effect and cannot be observed if the advection term of Eq. 1 is neglected. The advection patterns along which the particles move are relevant to fluid mixing.

The Stokes drift velocity field that represents these vortex patterns is a way to better visualize these traveling paths for greater time scale than the period of vibration. The vortex patterns were affected by wavelength and boundary conditions. For example, the right-most vortex was generated due to the rigid boundary along the right edge. The other vortices have sizes roughly comparable to the wavelength (4, 2 and 1 mm for panels C, D and E, respectively). The direction of advective flow was counterclockwise—along the wave propagation direction near the boundary and against the wave propagation direction in the middle. As fluid particles tend to circulate within vortices, long vortices are expected to enhance longitudinal mixing more than short ones. These velocity fields shown in Fig. 3C–E were used as the advection velocity (**u**_{D}) in solving Eq. 6.

### Advection and diffusion

To evaluate fluid mixing in terms of ion concentration, advection–diffusion equation (Eq. 6) was solved for different mechanical agitations. Diffusion in our study was nearly one-dimensional due to the geometry (i.e., the length dimension is dominant). The fluid homogenized in the radial direction within a fraction of second, while it took over an hour along the longitudinal direction. Considering this one-dimensional aspect, initial concentration was given so that it varied linearly along the length of fluid tube: (Cleft(x,t=0right)={C}_{L}+({C}_{R}-{C}_{L})x/L), where ({C}_{L}) and ({C}_{R}) are the initial concentration values at the left and right ends. Given values were 5 and 1 mM for ({C}_{L}) and ({C}_{R}), respectively. It should be noted that our results are based on normalized quantities. Thus, the absolute values of concentration do not affect our result, but the longitudinal pattern of concentration gradient matters. As the spatial grid for mass transport (advection–diffusion) was different from that for fluid dynamics, velocity values at grid points were interpolated from the fluid dynamics solution.

Over time, the fluid was mixed and homogenized due to diffusion and advection. To quantify the level of homogeneity, the mixing parameter *χ* was defined using the standard deviation of ion concentration across the fluid space:

$$chi left( t right) = 1 – {text{stdev}}left( {Cleft( t right)} right)/{text{stdev}}left( {Cleft( 0 right)} right).$$

(7)

The mixing parameter is zero initially and it approaches one as the fluid space becomes homogeneous. The absolute value of concentration does not affect the mixing parameter. Figure 4A shows the change of mixing parameter over time when the Corti fluid model is subjected to the mechanical agitation of Fig. 3C, and when there is no mechanical agitation (pure diffusion). When mechanically agitated, the fluid was homogenized quicker. To quantify the speed of homogenization, the mixing time was defined as the intercept between the initial tangent line to the curve of *χ*(*t*) and the line of *χ* = 1 (Fig. 4B). Initial slope was chosen to define the mixing time because the mixing at early stage characterizes the advective mixing well. In the mechanically agitated case, the mixing time was faster by up to one order of magnitude than pure diffusion (6.4 min for 80 m/s case versus 74 min for pure diffusion case, Fig. 4B). There existed a critical wave speed across which the dominant mode of mixing switched between diffusion and advection (Fig. 4C). For the considered micro tube and stimulation level, the critical wave speed was 7 m/s. The relationship between the mixing time and the wave speed is shown in Fig. 4C. The horizontal axis is also presented with the Péclet number (*Pe*), defined as (Pe={stackrel{-}{u}}_{D}h/D), where ({stackrel{-}{u}}_{D}) is root-mean-square Stokes drift velocity over the space, and *D* is the diffusion coefficient. The Péclet number represents the contribution in mass transport due to advection relative to diffusion. Consistent with the physical implication of *Pe*, the mixing time for advection–diffusion became inversely proportional to the wave speed when *Pe* > > 1. When *Pe* < < 1, the mixing time approached to the value of pure diffusion (filled square). Note that the physiological speed range (> a few m/s) corresponds to the domain of advective mixing.

### Physiological traveling waves and micro fluid mixing

Physiological traveling waves differ from the simple traveling wave of Eq. 3 in that the wavelength gets shorter as the wave propagates toward the apex (Fig. 5A). As a consequence, the wave speed (phase velocity) decreases as the wave propagates (Fig. 5B). The wave pattern of Fig. 5A approximates a measurement of chinchilla cochlea^{25}. The original measurement was made at a location where the best responding frequency was 9 kHz. Sound tones at different frequencies were delivered at 70 dB sound pressure level (SPL) and vibrations were measured at the basilar membrane (the bottom boundary of the OoC in Fig. 1A). From the frequency responses, the spatial vibration pattern was obtained using the scaling symmetry^{26}. A rough approximation we made was that the measured basilar membrane displacement is comparable to the cross-sectional radius change (*r–r*_{0}, in Eq. 3). Two recent observations may justify this approximation. First, there is significant phase difference between the top and bottom surfaces of the OoC (up to 180°^{11}) which results in cross-sectional area change. Second, the top of the OoC vibrates more than the bottom (basilar membrane)^{27,28}.

Using the same Corti fluid model as the one in Figs. 1, 2, 3 and 4, the stimulating boundary condition was replaced with this physiological traveling wave to obtain the results shown in Fig. 5. Stokes drift trajectories show fewer isolated vortices compared to the simple wave case in Fig. 3, despite multiple wave cycles over the simulated length. This resulted in fewer transport barriers that slow down longitudinal fluid mixing. The Stokes drift velocity sharply drops across the peak of wave envelope (x ≈ 3 mm, Fig. 5C). Instead of computing mixing parameter for the entire fluid domain, local mixing parameters in the tail (circle), peak (square), and head (rectangle) of the traveling wave were computed (Fig. 5D). The local mixing was considered for two reasons. First, we wondered in which region of traveling wave the fluid mixes better. As the wave parameters vary along the distance, the wave speed gets slower toward the apex, while the vibration amplitude peaks where the wave speed is slow (Fig. 5A,B). Second, cochlear fluid mixing must be a superposition of local events. Unlike the artificial case in Fig. 3, natural traveling waves do not fill the entire tube length—their wave envelopes have finite lengths. For instance, when the gerbil cochlea is subjected to a modest level of sound, 20 and 1 kHz sounds make approximately 2 and 9 mm-long wave envelopes, respectively^{29}. The envelope length depends not only on frequency but also on sound pressure—louder sounds result in wider wave envelope. Waves with longer envelopes have higher mixing parameters if the entire domain is considered. Local mixing matters because, in reality, Corti fluid mixing is an ensemble of many waves. At a specific location, the most efficient wave component will dominate the mixing. Therefore, mixing of the entire cochlear length by a single wave does not represent Corti fluid mixing well.

The fluid mixing was more active toward the tail (*x* = 0) of the traveling wave. This location-dependence of mixing is reasonable considering that the wave speed decays (from 100 m/s to less than 1 m/s, Fig. 5B), while the vibration amplitude does not change dramatically between the tail and peak region. The local mixing parameter does not increase monotonically in the head region of traveling wave (indicated with rectangle, Fig. 5C,D); the mixing parameter decreased over the first 10 min.