Oral cavity hydrodynamics and drag production in Balaenid whale suspension feeding

Balaenid whales feed on large aggregates of small and slow-moving prey (predominantly copepods) through a filtration process enabled by baleen. These whales exhibit continuous filtration, namely, with the mouth kept partially opened and the baleen exposed to oncoming prey-laden waters while fluking. The process is an example of crossflow filtration (CFF) in which most of the particulates (prey) are separated from the substrate (water) without ever coming into contact with the filtering surface (baleen). This paper discusses the simulation of baleen filtration hydrodynamics based on a type of hydraulic circuit modeling commonly used in microfluidics, but adapted to the much higher Reynolds number flows typical of whale hydrodynamics. This so-called Baleen Hydraulic Circuit (BHC) model uses as input the basic characteristics of the flows moving through a section of baleen observed in a previous flume study by the authors. The model has low-spatial resolution but incorporates the effects of fluid viscosity, which doubles or more a whale’s total body drag in comparison to non-feeding travel. Modeling viscous friction is crucial here since exposing the baleen system to the open ocean ends up tripling a whale’s total wetted surface area. Among other findings, the BHC shows how CFF is enhanced by a large filtration surface and hence large body size; how it is carried out via the establishment of rapid anteroposterior flows transporting most of the prey-water slurry towards the oropharyngeal wall; how slower intra-baleen flows manage to transfer most of the substrate out of the mouth, all the while contributing only a fraction to overall oral cavity drag; and how these anteroposterior and intra-baleen flows lose speed as they approach the oropharyngeal wall.


Introduction
All members of the Suborder Mysticeti use baleen to separate the prey (particulates) from the sea water (substrate). However, the manners in which the filtration process takes place turn out to be quite different in each family, leading to different ways of generating drag. Balaenid whales, including the bowhead whale (Balaena mysticetus Linnaeus 1758) and right whales (Eubalaena spp. Linnaeus 1758), represent one of few groups of large-bodied organisms that a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 use crossflow filtration (CFF) as a means to collect large amounts of small and slow-moving prey [1][2][3][4][5][6][7][8][9]. Baleen disposed longitudinally within the oral cavity is used to separate prey from the water, which flows continuously in and out of the mouth and filter due to the whales' forward motion (Fig 1). The filtration is deemed as "crossflow" because most of the particulate moves parallel to the filtering surfaces (baleen) rather than into them as with "throughput" filtration [9,10]. Here the prey directly accumulates at the oropharynx rather than being brought there by the tongue after sieving; in other words, a balaenid whale feeds on a suspension rather than a sieved mass tangled in the baleen and fringes after skimming [9]. This is crucial here, since the prey (primarily copepods) turn out to be significantly smaller in size than the width of the intra-baleen gaps forming the filtration surface [6,10].
In contrast, and again during filtration, the baleen of rorquals (Balaenopteridae) are exposed not to the open ocean (except momentarily, during engulfment) but to the intraoral water-prey mixture engulfed during a lunge and contained within the closed buccal cavity [11][12][13]. And so the movement of water past and out through baleen is instead enabled by internal pressure buildup within the buccal cavity, as caused by the contraction of musculature embedded within the ventral groove blubber (VGB) [14,15]. This intermittent filtration mechanism is also used by the gray whale, Eschrichtius robustus, although in gray whales the engulfment of a single mouthful of prey-laden water occurs via intraoral suction generation rather than by rapid lunges, and there are far fewer throat grooves to expand and hold engulfed water; relative to rorquals there is less intraoral pressure generation and less gular musculature to contract and expel the excurrent water past the baleen filter [1,7,16].
A quantitative description of these intra-oral flows for any of these mysticetes has remained elusive, however, due to the obvious impossibility of studying such large animals in the laboratory. Only by parsing the problem into parts amenable to experimental and theoretical investigation, and along with using kinematic data collected in the field, can we hope to further understand what is bound to remain out of view. Moreover, exhibiting what appears to be the least dynamic of the filtration strategies used by baleen whales, balaenid CFF offers the simplest flows available for study.
As a first step towards this goal, a prequel paper [9] has discussed results of a flow tank study of baleen hydrodynamics near and through a section of real baleen, showing CFF to arise mainly from their cambered shape. This process is also enabled by the anteroposterior and mediolateral pressure gradients that split the current entering the mouth into a lateral intrabaleen (IB) flow and an anteroposterior (AP) flow moving towards the oropharyngeal wall (Fig 1). In the flow tank study the AP pressure gradient was created externally by the flume's pump system, thereby recreating in a very approximate manner the lower pressure found near a whale's Posterior Opening (PO). The AP flow is the faster of the two flows, thus allowing the transport of most of the engulfed slurry of prey towards the oropharynx. IB flow, whose sole purpose is to reduce the water contents of the slurry, is significantly slower but capable (simultaneously) of transferring a large mass of water, owing to the large surface area of the baleen racks in comparison with that of the mouth inlet (Figs 1 and 2; see also Fig 1 and 2c in [9]). The great efficiency of balaenid CFF stems from the generation of slow IB flows in conjunction with fast AP flows, which minimizes the number of prey items that pass through the filter or get entrapped in the baleen plates and fringes. Moreover, slow IB flows ensure that the time needed to clog the baleen and fringes (when happening) is much longer, in comparison with conventional "dead-end" sieving scenarios [9,10]. Finally, the study shows the robustness of CFF in taking place regardless of the precise emplacement of the tongue and lip walls.
Flow tank studies of balaenid baleen systems are by necessity limited in scope given the truly large sizes of these animals. With actual oral cavities measuring up to 3-5m in length  (Table 4); A IBchannel , Surface area of each intra-baleen canal (or channel) ( Table 4); A PO , Surface area of the Posterior Opening (Fig 2) ( Table 4) Total drag force; h HT , Average height of the baleen plates along each rack (Fig 2); k canal , Friction coefficient; k APL , Friction coefficient for APL canal segments; k APT , Friction coefficient for APT canal segments; k IB , Friction coefficient for an IB canal; L body , Body length; L mouth , Length of the oral apparatus; equal to N b d baleen ; M c , Body mass; N b , Number of baleen plates on one rack; P ext , Pressure outside and near the Posterior Opening; P i , Pressure in the APT canal at site i ( Fig 3); P in , Pressure at the mouth's entrance; P pi , Pressure in the APL canal at site i ( Fig 3); P PO , Pressure at the Posterior Opening (PO); Re, Reynolds number (= speed x length/ν); S wetted , Wetted surface area of the body (closed mouth); U ext , Speed of seawater flow external to the body and near the PO; U in , Flow speed entering the mouth (sub-rostral); U i , APT flow speed at station i; U pi , Inlet intra-baleen flow speed; at station i; U Li , APL flow speed at station i; U whale , Travel speed of a whale during foraging; x IB , IB friction coefficient parameter (Eq 13); y rostr , IB friction coefficient parameter (Eq 13); z rostr , IB friction coefficient parameter (Eq 13); δ(x), boundary layer thickness a distance x along canal wall; δ IB , boundary layer thickness at the outlet of an IB canal; δ APT , boundary layer thickness at the outlet of an APT canal segment; γ, Drag correction factor due to depth; ν, Kinematic viscosity of sea water; ρ w , Density of seawater; AP, Anteroposterior; APT, Anteroposterior, lingual side; APL, Anteroposterior, labial side; BHC, Baleen Hydraulic Circuit; CFF, Cross-flow filtration; IB, Intra-baleen; PO, Posterior Opening.
(anteroposterior), 1-3m in height (dorsoventral), and 1-3m in width (mediolateral), the typical size of flumes (3m 3 or less) makes hydrodynamic analysis of a fully-sized baleen rack all but impossible. Thus only small sections of baleen can be studied experimentally, to yield information with regards to the local flows near and through baleen [5,9]. However, and thanks to a rack's long serial axial symmetry where pressure and (flow) celerity gradients vary little over length scales of 0.3m or less (Fig 1), and where the camber and triangular shape of individual plates also vary little along the length of the whole baleen rack, such data can be used in mathematical modeling of the global flows characterizing the entire oral cavity. This has been done here via a hydraulic circuit approach commonly seen in microfluidics [17], but modified for the high Reynolds number flows relevant to balaenid hydrodynamics. A so-called Baleen Hydraulic Circuit (BHC) model was devised to calculate the pressure drops and flow speeds  figure), blue arrows indicate direction of water flow though and around baleen filtering apparatus in life as well as in experimental flow tank trials and computational modeling calculations (hypothetical but predicted from data of current study and previously published experiments [3,5]). Water can flow anteroposteriorly (AP) within mouth along the tongue (APT channel) or the lip (APL channel). Filtered water exits the mouth via paired posterior openings (PO). Oropharyngeal opening which leads to esophagus lies near oral floor caudal to the tongue root. Dashed red box indicate the approximate location of the shortened mini-rack studied in a previous flow tank study [9]; dotted green box shows the system under consideration in this paper.
resulting from the viscous friction present in between the intra-baleen gaps and the body drag that ensues. This is a low-spatial resolution model, but simple enough to be solved by common spreadsheet programs. Its inputs includes morphological data characterizing the dimensions of the various canals defining the oral cavity (Figs 2 and 3), flume data to set the proportion of fluid moving into a particular baleen canal in relation to the mass that continues past it [9], and modeling of the viscous friction between the flows and oral tissue. These BHC calculations first yield an estimate of the extra drag force arising from baleen exposure to the open ocean, which in turn is used to infer oral cavity hydrodynamics. As shall be seen here, and in comparison to that of non-feeding travel, baleen exposure can triple the net wetted surface area of a whale, and increase total drag by 100% to 400% depending on mouth inlet area, or, in other words, on how much a whale lowers it mandibles and/or cant its lips (Fig 2). These calculations turn out to agree reasonably well with the drag inferred from the kinematic data collected with digital tags deployed on feeding whales in a glide [18].
Ultimately, knowing how much drag there is adds new insights in the trade-off between the energy gained through ingestion of prey and the energy lost through the metabolic expenditure of foraging. All of this has ecological implications, particularly with regard to the minimal copepod densities (~400/m 3 ) and the not-to-exceed swim speed (~2m/s), needed to assure surplus harvest yields which on a yearly basis must guarantee survival over long fasting periods. These aspects and more shall be addressed in a sequel paper.

Materials and methods
The physical problem Efficient filter design begins with packaging a large filtering surface area within an enclosure that generates the least drag [6,9]. This is no different with balaenids, as exemplified by the case of a 10m, 15,155kg whale simulated here (Tables 1-4), for which longitudinally-oriented oral tissue exposes 144m 2 of surface area to the open ocean. This is a large area in comparison calculated in the BHC model. Arrows represent the cross section-averaged flow speeds (the longer the arrow the faster). On the lingual side, anteroposterior flow decelerates because of the mass loss incurred by intrabaleen flow. On the labial side, flow accelerates from the mass gained from intra-baleen flows. Variables A in and A po correspond to surface areas of the gap between the baleen rack and tongue's side wall and to the posterior opening (Fig 2).
https://doi.org/10.1371/journal.pone.0175220.g003   tongue and oral cavity floor surface, albeit imperfectly). With such a large exposed filter comes the viscous friction generated by the boundary layers attached to these surfaces. This, in turn, significantly increases a whale's drag, and to levels exceeding twice the drag sustained during non-feeding travel [18]. The flow dynamics within the oral cavity takes place over a wide range of scales. First, fluids need to move over the entire length, height and width of the mouth, which span scales of one meter up to several meters as discussed previously. Water travels similar distances over the entire length ("span") of baleen (1 to 2m), but also over much smaller distances along the few centimeters characterizing the width (or "chord") of baleen (7cm), and through the gaps separating the baleen (1cm) [13]. With foraging swim speeds spanning a range of 0.3m/s to 1.2m/s [18,[19][20][21][22], and the generally decreasing internal longitudinal flow speeds towards the oropharynx [9], the Reynolds numbers associated with the latter stand in the range of R e~1 0 5 -10 6 ( Table 3). These contrast with intra-baleen flows where R e~4 x10 3 to R e~1 x10 4 owing to the much smaller pipe diameters involved (Table 3).
(These values are at least twice as large as those characterizing the Hagen-Poiseuille pipe flows of many biological systems [23]). Yet smaller scales of 1mm or less are found in the boundary layer and turbulence present over and near the baleen plates, as well as over both tongue and lip walls.
Studying in the laboratory and numerically the effects of length scales spanning up to five orders of magnitude is difficult. Generally, Computational Fluid Dynamics (CFD) techniques are powerful and accurate enough to handle such a range. But these would require unaffordable costs in computer hardware and time with the problem at hand, since tripling the surface area exposed to flow would, in comparison to the already-challenging simulations of closedmouth swimming [18,24], require approximately five times more grid points. (In typical CFD simulations, the bulk of the computational domain's grid points are found near and within the fluid's boundary layer coating the exposed surfaces. Thus tripling the latter would in effect increase the associated volume by a factor~area 3/2~33/2 ). This is problematic given the substantially large input parameter space explored here (Tables 1 and 2). We note that attempts at using CFD to model balaenid filtration and drag have been described in [18]. Our approach has been to develop a simpler hydrodynamic simulation scheme, but one that is complete and fast enough to provide a first look into the relevant hydrodynamics, thereby providing a road map for future CFD and flow tank studies.
The BHC is a low spatial-resolution hydraulic circuit model capable of calculating the local pressures and flow speeds both anteroposteriorly and mediolaterally and through the baleen/ fringe complex, while accounting for most of the viscous friction sustained by these flows. Note that using hydraulic circuit modeling excludes the existence of large scale vortical motions taking place over the largest scales of the mouth (i.e.,~1m), as has been invoked in other organisms [25]. Conceptually, this means that overall flow celerity will be driven by pressure gradients rather than inertial effects, as applied over the entire length, width and height of the mouth apparatus. Herein the main anteroposterior pressure gradient arises from the low pressure at the PO associated with the external flows accelerating around the oblong geometry of the body [3,9]; and the medio-lateral pressure gradient arising across the racks due to baleen curvature [9]. The BHC model is described qualitatively in the following sections, with the detailed mathematics applied to the simpler and pedagogical 4-baleen system discussed in the six sections that follow the Final Remarks (Modeling Details [1][2][3][4][5][6]. Application to the more general N b~3 00 baleen system is discussed in Modeling Details 7, which is followed by a list of the mathematical symbols and acronyms.

Simulated oral apparatus
BHC modeling calculates the flows and pressures at the locations shown in Fig 3, including within each intra-baleen (IB) gaps, as well as in the anteroposterior gaps/canals, i.e., both lingual (APT) and labial (APL). Several sites are used for the calculation of both flow speed and pressure, while others are assigned to either quantity. The pressures P i and P pi (at site i along a rack) correspond to the pressures in the APT and APL canals respectively, and by extension to the pressure at both ends of the IB channel separating both AP canals there. Obviously, the rectangular cavity of Fig 3 is a rough approximation, in comparison to the actual, curved 3-dimensional cavity exemplified in Figs 1 and 2. This applies to the baleen system as well, where racks of straight, equal-span baleen is simulated as an approximation of the actual racks of curved and unequally-spanned baleen [16]. Although leading to a drastically simplified description of the fluid dynamics, the model cavity is based on the averaged dimensions found in vivo, and should yield useful estimates without losing too many of the important details. Note that in contrast to the simple system of Fig 3, the tongue of a whale does not present a boundary to the flow over the entire height of the mouth (Fig 2). Whether this is a good approach will depend on future studies of tongue emplacement and volume during feeding.
Another approximation involves calculating the viscous friction associated with an APL flow that origins entirely from the lateral motions (through baleen) of fluid initially located in the APT canal. This calculation omits altogether the contribution of the water that has entered the APL canal through the anterior opening of the lips (compare Figs 1 and 3). The effect that such additional flow has on overall APL flow and total body drag is currently unknown, and probably assessable only with more detailed CFD simulations. At present, it is suspected to be small, again due to the small of amount of extra fluid mass such opening allows in, compared to the mass flowing through the APT canal.
Finally, and given the assumed rectangular geometry of a BHC oral cavity, the resulting BHC solutions are equivalent to 3-dimensional flow moving horizontally in the bulk. In other words, the flow is the same dorsoventrally, in contrast to what may happen near a real section of baleen [9].

BHC flow speed laws
The model enforces conservation of the flow mass rate entering and emerging from each canal (Fig 4) [9,26]. It follows that the water entering the mouth and moving anteroposteriorly loses mass to the intra-baleen flow (Fig 1) and progressively decelerate as it approaches the oropharyngeal wall [9]. Such speed reduction occurs regardless of the existence of the pressure gradient arising from the high pressure at the mouth inlet and the low pressure at the posterior opening [3]. (In a single-outlet pipe, such a gradient would accelerate the flow). This trend of slowing APT flows has biological implications, including a slow-moving food slurry that is easier to swallow.
Calculation of the flow speeds are obtained after solving a system of equations implied by mass flow rate conservation. Details are shown in Tables 5 and 6, and Eqs 26-28 in Modeling Details 7. But complete solution isn't achievable until a second set of equations connecting APT and IB flow speeds is defined. Used here is the flow-splitting coefficient C, which calculates the ratio of speeds within an intra-baleen canal (U IB ) to the APT speed near the inlet of the given IB canal (U AP ), or C = U IB /U AP . This coefficient can be measured in flume studies and is considered as constant in both time and space (i.e., over all IB canals). Values were extracted from data of the prequel paper [9], namely C~0.44 ± 0.04 and 0.53 ± 0.03, obtained with and without the prey particles in the rack respectively. Interestingly, these results are independent of tongue wall locations (or APT canal width D in ) and over the entire range of the swimming  Tables 1 and 2. https://doi.org/10.1371/journal.pone.0175220.g004 Table 5. Flow rate conservation (left) and flow splitting (right) equations in the 4-baleen BHC system. Indices correspond to the locations shown in Fig 3.

APT canal sections IB canals
https://doi.org/10.1371/journal.pone.0175220.t005 Balaenid oral cavity hydrodynamics and drag speeds used by foraging balaenids [14][15][16][17], an effect hypothesized to arise mostly from baleen camber [9]. Untested, however, is the dependence of C on the intra baleen separation (d baleen ), thus limiting our sensitivity study to a narrow range (d baleen = 0.01-0.02m). The ratio U IB /U AP leads to the equations in the right column of Table 5, which after merging with those of the left column, yield the speed solutions shown in Table 6.

BHC pressure drop laws; canal friction and boundary layer
The pressure differences at the ends of each canal are estimated with a type of pressure drop equations commonly used in constant-diameter pipe flow engineering [26], adapted here to the rectangular canal geometry implied in Fig 3. Being sustained in pipes/canals characterized by Reynolds numbers well-exceeding Re~2000, the pressure gradient (ΔP) across both ends of a given canal is connected to the square of the flow speed at its inlet, i.e., ΔP = k canal ½ ρ U inlet 2 . Sometimes known as the Darcy-Weisbach equation, this relationship differs from the wellknown Hagen-Poiseuille law where ΔP is only proportional to U inlet [23]. Application of this formula is shown in Table 7 and this is where the BHC accounts for the viscosity and turbulence effects affecting the drag. Parts of the canal friction coefficient k canal , which varies at each canal, is calculated via the modeling of the flows as undeveloped flows, where the maximum lateral extent of the boundary layer coating the canal (flat) walls remain smaller than a canal half-width [26]. Moreover, and as further explained in Modeling Details 4, the coefficient is further adjusted to insure the absence of adverse pressure gradients within the cavity, along with pressure levels at the PO that are consistent with those exterior to the body.

Table 6. APT and IB flow solutions of the 4-baleen BHC system. Indices correspond to locations in Fig 3.
A IBchannel is the cross section area of the intra-baleen space and equal to the product h Ht c baleen .

Viscous friction drag of the mouth, versus total drag
The drag produced by the oral cavity and body is an important physical characteristic that provides the means of estimating the metabolic expenditures connected with swimming and feeding [27]. But calculating the drag from BHC output also allows for checking its accuracy since total body drag can be inferred from digital tag data collected in the field [18]. Drag measures the resistance encountered by a body moving through a fluid [23,26]. In steady motions the common sources of drag are friction drag and pressure drag, with the latter arising from the turbulence of the whale's near-wake [28], and the former from the viscous friction experienced by the fluid moving past the solid boundaries of exposed skin and baleen. Total drag is thus calculated from the sum of the mouth viscous friction drag (F D mouth ) associated with the flows passing through the oral cavity, and the friction and pressure drag combination created by the (external) flows moving around the body (F D body ): where and Eq 1 is an approximation that neglects the interaction between the external flow and the internal flow emerging out of the PO. Eq 2 corresponds to the standard drag equation used to model the effects of water resistance during non-feeding travel [18,[29][30][31][32]. The factor in brackets is an approximation of the wetted area of the body, minus that of the oral cavity inlet area. Even with an open mouth, most of a whale's body-including the lips-deflects fluid around its exterior (Fig 1). The only area to subtract is the body feature that does not deflect fluid, namely that of the oral cavity inlet (2A in ) (Fig 2). Parameters C D closed mouth and S wetted correspond to the non-feeding travel drag coefficient and wetted area of the body respectively, both in closed mouth configuration. The former is set at 0.0035 (15m body length) and 0.0059 (8 and 10m) per the glide study of Nousek McGregor [18], and the latter calculated using a correlation discussed in [29] The factor M c corresponds to a whale's body mass (kg) and is estimated from a correlation established in [33]: Finally, the factor γ represents a swimming depth correction factor that takes into account the (body) drag increases incurred when swimming near the surface [29,34]. Herein γ = 1, corresponding to foraging bouts taking place at depths exceeding three times a whale's body's maximum (vertical) width (as in [18]).
Eq 3 describes the mouth friction drag term arising from the viscous and turbulent energy dissipated within the oral cavity. The first term is proportional to the total energy density-i.e., kinetic (½ ρ w U 2 ) plus potential (pressure)-of the flows entering the mouth through area A in at speed U in , minus that of the flows exiting the mouth through the PO at speed U PO . The difference measures the dissipative losses, viscous friction and turbulent, incurred by the flows while moving in the anterio-posterior direction (see Modeling Details 5 for details). The second term represents the energy dissipated when the APT flows are redirected into IB flows, moving through baleen, and then redirected into APL flows. These effects are modeled here by looking at a rack of baleen as one of closely-spaced, high-angle-of-attack hydrofoils of known drag characteristics (C D baleen = 0.25 here). This second term is calculated during input data entry and is used in the final adjustment of the IB canal friction coefficients k IB . A detailed derivation is described in Modeling Details 5.
The last ingredients are the drag coefficients C D total and C D mouth corresponding to total body drag and mouth-generated drag obtained from Eqs 1 and 3 respectively: Usually the wetted area of the entire oral cavity would be the preferred reference surface area for the calculation of C D mouth . However, using the same (closed-mouth) body reference surface area S wetted for both drag coefficients will permit their direct comparison when assessing the importance of mouth drag in comparison to that the drag generated by the rest of the body.

BHC input parameters and constraints
Input parameters and constraints for BHC model calculation involve most known factors relevant to balaenid suspension feeding (see Figs 2 and 3 and Tables 1 and 2): • Basic fluid dynamics: sea water density (ρ w ), kinematic viscosity (ν), ocean static pressure measured at the surface (P atm ) and whale swimming speed (U whale ); • Morphology: mean baleen chord (c baleen , assuming consistent serial axial morphology of plate geometry and camber, as mentioned above), spacing along rack (d baleen ), and total number of plates on the rack (N-1); • Width (D in ) of the APT canal between tongue wall and medial side of a baleen rack; • Lateral width of posterior opening (w PO ) characterizing spacing of the APL canal between lip and labial end of the rack (See Eq 9 below); • The baleen plate span (h HT ), averaged over the length of a rack, used for filtration during foraging; • Flow splitting coefficient (C), as extracted from flume data [9]; • Coefficients f APT , x IB , y rostr and z rostr , used in the friction coefficient profiles calculated via Eqs 13-16 (see Modeling Details 4). Their values are adjusted to: 1) produce internal pressures at the PO that are consistent with pressures immediately external to it; 2) yield internal pressure profiles free of adverse gradients; and 3) insure viscous energy dissipation in between the baleen plates similar to that of the hydrofoil model used in Eq 3. This adjustment is described in Modeling Details 1 and 4.
Parameters D in , w PO and h HT characterize morphological attributes that can be evaluated from field observations (ranges stated in Table 2), and likely to vary with the size and density of a prey patch. The calculations that take place during input data entry are further explained in Modeling Details 1.  Table 5) which are typically high, but erroneous given the lack of knowledge on the size and geometry of that area. The figure shows both APT and IB flows becoming slower further down the rack, a direct result of the APT flows losing mass to the IB canals per conservation of the mass rate [9]. Greatly reduced speeds near the oropharynx are desirable as they will facilitate the swallowing of the food slurry accumulating in that area. The labial AP flows (APL) are shown as well, and are calculated from Eqs 24 and 25 (Modeling Details 7). As expected, this flow increases in speed as it moves along the canal, after picking up the fluid mass exiting the intra-baleen channels. The figure shows a similar comparison with the pressure drops sustained. Along with expectations, the AP pressures are greater than the IB pressures over the entire length of the rack. Moreover, both monotonically decrease to similar values near the oropharyngeal wall. This pattern thus excludes backflows in any one of the canals located upstream, as they are inconsistent with the flow solutions of Table 6 and Eqs 27 and 28, and unlikely to occur in a real whale. Finally, each calculated pressure profile is continuous within numerical discretization variations, to yield the same pressure difference between any two points in the cavity (within a fraction of a percent), as calculated from integrating pressure increments along arbitrary paths in the APT and/or APL canals starting and ending at those same two points.

Pressures and flow speeds
Typical values of the friction coefficients k canal are shown in Fig 5 (calculated from Eqs 13-16). Generally, the coefficients are large where the boundary layer is "thick" (relative to pipe diameter) and/or when the flow speed is small. The values obtained for the IB flows are quite similar in order of magnitude to those encountered in microfluidic and hydraulic flows moving through bends of various angles [35,36]. The friction coefficients of the wider APT and APL canals are significantly smaller in comparison but also comparable to the values encountered in similarly-sized straight-pipe and nozzle flows of similar Reynolds numbers [26,36,37].
Flow speeds obtained at different values of the flow-splitting coefficient C are shown in Fig  6. Different values of C correspond to various states of baleen canal obstruction by the prey tangled with the baleen fringe fibers [9]. Clean baleen would be characterized by C~0.50, and C~0.40 by obstruction levels likely to occur with particulate densities found in the field. Lower values are possible and presumably correspond to yet higher levels of rack obstruction.
The figure shows that the lower the C, the slower the IB flow in the anterior third section of the baleen rack, but faster in the posterior two-thirds. This pattern is necessary to keep up with the system's in-flow mass rate which is held at the same value in the figure. In the APT canal, lowering C translates into faster speeds along the entire rack. To put it in another way, both flows would be concentrated along the anterior portion of the baleen rack at the higher values of C, but become more evenly distributed along the rack at lower values. Such a pattern suggests that baleen clogging would occur at first along the anterior section of the rack, to progressively expand further down the rack over time. This isn't happening here, with the BHC being essentially a steady-state simulation of the flows. Another important aspect that has been left out of the figure is the decreasing mouth inlet flow rate U in and higher inlet pressure P in that would also accompany increased baleen clogging. Implementing this scenario over time would have resulted in the gradual lowering of the APT and IB flows of What the simulation shows, however, is that tongue emplacement is the more important factor, not only for flow speed, but also for drag generation as shown next.

Mouth viscous friction drag
Values of the mouth friction drag calculated via Eq 3 are displayed in Figs 9 and 10, both versus the square of the swim speed. The former is a comparison among the three body sizes investigated, namely 8m (with a 2.4m-long oral cavity), 10m (3m-long) and 15m (4.5m-long). All three sizes feature different baleen spacing in the range of 0.01-0.02m (Table 1). Fig 10 shows the mouth drag generated by a 10m case in which three APT canal widths have been simulated (D in = 0.30, 0.50 and 0.72m). All simulations were performed with 299 baleen plates and at value of the flow-splitting coefficient C equal to 0.44 ("with prey" [9]).
Comparing these figures show an unambiguous proportionality with U whale 2 for all cases investigated, with the slopes strongly dependent on the value of D in in Fig 10, and on body size in Fig 9. With respect to the latter, Eq 3 suggests that these variations in slope have more to do with the value of the inlet area D in h HT than with body size. Here the values of D in h HT turned "with prey" [9]) and C = 0.53 (dotted; "no prey" [9]). All three cases involve the same mouth inlet flow rate. The rest of the inputs used in both frames where: 10m body length, U whale = 1.0m/s, d baleen = 0.01m, c baleen = 0.07m, h HT = 1.5m and D in = 0.5m. Further details are listed in Tables 1 and 2. https://doi.org/10.1371/journal.pone.0175220.g006 out at 0.48m 2 (8m body length), 0.75m 2 (10m), and 1.69 and 1.54m 2 (15m; case A and B respectively). In Fig 9 (10m body length) these inlet areas varied in a similar fashion, namely, at D in h HT = 0.45, 0.75 and 1.08m 2 . These results show that inlet cross section area is the dominant factor in oral cavity-generated drag across body sizes, with the implication that drag-control during foraging is to be mainly achieved by mandible and lip movement. The slopes listed in the captions of Figs 9 and 10 should allow the interested reader to compute mouth drag at various speeds without ever running the BHC. But another, and more general correlation can be derived from the numerical data and Eq 3 as follows. Fig 11 shows the ratio of the first term "T 1 " within the curly brackets of Eq 3, here normalized as 2T 1  is associated with the variations that occur while tuning the f APT coefficient). Using this value in Eq 3 permits complete reconstruction of the BHC-calculated mouth drag, i.e., with racks of 300 baleen each, and with the morphology not too far removed from Table 1. More importantly, it also permits the understanding of the scaling of mouth drag with respect to all of the relevant morphological inputs: for example, of T 1 scaling linearly with D in , in contrast to the second term of Eq 3 ("T 2 ") scaling with the cube of D in ; or with T 1 being largely insensitive to baleen spacing (d baleen ), versus T 2 which is inversely proportional to d baleen 3 (via L mouth defined as L mouth = N b d baleen ).
Most interestingly, such differential scaling ensures that about 60 to 90% of the mouth drag is generated by the T 1 -term, i.e., AP flows energy dissipation, with the T 2 -term representing that of the baleen and being important only in cases of large D in (> 1m) used in conjunction with small d baleen (< 0.01m). (In such cases T 2~5 0% of the total mouth drag). Even though baleen exposed to seawater triples the wetted area of the body, it appears that their contribution to the drag turn out to be more modest. Per the BHC, this would come, mainly, from the   Another scaling property of Eq 3 is the mouth drag increasing linearly with the average length of baleen (h HT ). This scaling is trivial here since both AP and IB flows are defined to move strictly in the horizontal (coronal) plane of a rectangular cuboidal cavity of uniform depth. On the other hand, scaling with respect to flow-splitting (C) appears minor in comparison, yielding 10% changes in mouth drag in the range C~0.3-0. 5.
In comparison to the 100-165N total drag sustained during non-feeding travel at 1m/s [12], the viscous friction of the mouth turns out to be as large if not larger, as with the case of the 8m (F D mouth 130N), 10m (211N) and 15m whales investigated here (444N with case A; 435N with case B). Thus mouth drag and mouth drag coefficient (C D mouth ) end up doubling, tripling and even quadrupling total drag, depending of the amount of mouth inlet area. This trend is shown in Fig 12 which the coefficient of total drag C D total (calculated from Eq 6a) versus swim speed. This data is also compared with those of a field study of foraging and traveling right whales [12]. Both sets of data (numerical vs. experimental) appear in agreement with respect to the drag-range that can be generated, presumably through the variation of mouth inlet area (D in h HT ).

Energy dissipation pattern in the IB canals along the baleen rack
A rack of baleen can be seen not only as an assemblage of hydrofoils, but also as a parallel network of canals each independently dissipating viscous energy. From the latter perspective emerges an alternate approach to estimating energy dissipation and attendant drag codified in Eq 20 below. As further discussed in Modeling Details 5, each canal (located in between baleen plates i and i+1) dissipates viscous energy at a rate given by D in h HT U pi x k IB (i) x ½ ρ w U pi 2 . Such a rate (at each baleen station) is shown in Fig 13 and turns out to be maximal along the anterior ) versus swim speed and body size. Calculated from Eqs 1-3 and 6a and for C = 0.44 ("with prey" [9]). Shown are the following body sizes: 8m (blue triangles); 10m (all black symbols, with D in = 0.50m (crosses), D in = 0.72m (times) and D in = 0.30m (starbursts); and 15m (yellow diamonds for case A and brown diamonds for case B). This numerical data is compared with the drag coefficient calculated from kinematic data collected during feeding (red filled circles labeled "F") and nonfeeding transport (red filled circles labeled "T") [18]. Note that the foraging swim speed of balaenids has been observed in the range of 0.6 to 2m/s [18,[19][20][21][22].
https://doi.org/10.1371/journal.pone.0175220.g012 third of each rack. Although more difficult to quantify, the canal sections that forms the APT canal would show an analogous pattern, if only for the reason that APL flows are the fastest there as well (Fig 4).

Prey collection and drag control
The hydrodynamic modeling described here has highlighted two novel aspects of balaenid cross-flow filtration hydrodynamics, namely the pattern of decelerating APT and IB flows towards the oropharyngeal wall; and the substantial increase of body drag, as contributed by the mouth's viscous friction by flows going past the lingual and labial walls as well as through a very large number of tall and narrow baleen canals. The connection between these is key to understanding how the collection rate of prey, as determined by the rate of seawater intake D in h HT U whale , can be balanced against the generation of drag and attendant metabolic expenditures in the energy management of foraging. Thus in-flow rate is adjusted by the separate or combined use of the swim speed and mouth opening area, the latter achieved via the lowering of the mandibles along with the canting of the lip walls (Fig 2). All three parameters are also major determinant of drag (Eq 3 and Figs 9, 10 and 12), but in different proportions. Doubling the swim speed to collect twice as many prey items (per unit time) would lead to quadrupling the drag and required fluking thrust. As further explored in a sequel paper, there will be a speed "of diminishing returns" beyond which the collected prey will not yield enough energy to compensate for the metabolic activity of harvest. Interestingly, BHC modeling suggests that increasing mouth area instead of increasing the swim speed is the better strategy, as both drag and in-flow rate are both proportional to mouth area. The narrow range of swim speeds observed in the field, namely 0.6-1.2m/s [18,[19][20][21][22] would certainly be consistent with this hypothesis.

Comparative biomechanics of drag production
Balaenids and rorquals differ in how they create drag during filtration. With the former, drag is caused mainly by the viscous friction sustained by water moving anterioposteriorly past the lingual and labial walls, and to a smaller extent, through the baleen canals. Exposing about 600 baleen plates to ocean flows ends up tripling the body wetted area and double or triple the drag sustained with a closed mouth (Fig 12) (See also [18]). With rorquals, body drag during filtration is caused entirely by the friction and pressure drag generated over the tadpole-like body shape that results from the distension of the filled buccal cavity. What is different here is the water emerging out of the distended cavity originating from within the mouth rather than from the open ocean, thus contributing very little to the external flows that move about the body (to generate drag). Even in the shape of a tadpole, and while gliding following a lunge, the level of body streamlining is high enough to lead to comparatively lower drag coefficients, i.e., in the range of C D~0 .003-0.005, as with other streamlined animals or solid bodies [30,37]. Interestingly, rorqual body drag during filtration is quite small in comparison to that sustained during a lunge, typically smaller by a factor of 3 [27]. In other words, when taking into account the drag sustained over foraging in addition to filtration, both balaenids and rorquals sustain three-fold increases in drag, in comparison to that of non-feeding travel.

Large body size is always better for CFF
Increasing cross flow filtering efficiency begins with delaying the onset of filter clogging by reducing the flow speeds across the filtering surface. Speed reduction can be obtained by lowering the lateral pressure gradient, or better, by increasing the area of the filter in relation to that of the inlet [9]. The latter follows not only from BHC modeling, but more generally from mass flow rate conservation within the mouth and baleen. Further, and as mentioned in Results (sub-section Mouth viscous friction drag), slow intra-baleen flows reduce the drag associated with baleen, thereby enhancing the energy economy of large filtering surfaces.
Another (and less obvious) way of increasing efficiency is to reduce sea water contents of the prey/water slurry approaching the oropharyngeal wall, in order to lower the osmotic load of the kidneys, as well as to reduce speed of the nearby APT flows to facilitate swallowing a bolus of accumulated prey. Both feats are accomplished simultaneously since slower APT speeds occur at the end of both racks, again by virtue of mass rate conservation and by having a very large filtering surface in comparison to the mouth inlet area (Fig 4). (The filtering surface areas is 9m 2 in comparison to mouth area of 1.5m 2 . Here A filter = 2 h Ht 300 d baleen , with h Ht = 1.5m and d baleen = 0.01m for a 10m whale using 300 baleen plates per rack; see Table 1 and Fig 2; and A mouth = 2 D in h Ht , with D in = 0.5m). As further discussed in Modeling Details 6 (Eq 21), the ratio a of APT flow speeds at the last baleen plate, to that at the rostrum (or a = u 300 /U whale in a 300 baleen-per-rack system) can in fact be estimated by a~(1-CA IBfilter /A in ) Nb (with N b = 300). Clearly, the latter (and Fig 14) implies slower oropharyngeal flows with larger filtration area, which in turn supports the idea of large body sizes promoting greater CFF efficiency.

Final remarks
The BHC model is a simplified description of oral cavity hydrodynamics in balaenid whales. Albeit of low-spatial resolution, this model takes into account the viscous friction at work, particularly through the baleen canals which substantially increases a whale's wetted surface area. The model relies on several measureable morphometric inputs, along with hydrodynamic data obtained from a flume study of baleen rack [9]. BHC modeling yields a mapping of the pressures and flow speeds, both mediolaterally and anteroposteriorly, and an estimate of the drag generated by the mouth's viscous friction. So far, the calculated total drag agrees well enough with tag data [18]. Finally, the modeling has highlighted the role of large body size in enhancing filtration efficiency.
With the essential links between body drag and the patterns of intra-oral patterns flows now identified, further study is needed to clarify the role of several factors unmentioned here, preferably using Computational Fluid Dynamics (CFD) simulations of the model cavity. Most useful would be the clarification of the role of the flows entering the two gaps at the mouth's inlet dorsal to each baleen rack (currently closed), as created by the canted lips and lowered mandibles (compare Figs 2 and 3). The results discussed herein suggest their overall effects on drag as minor; but these could be important enough for the correct determination of the friction coefficients of the IB canals along the anterior portion of the racks, perhaps eliminating altogether the need for their adjustment ("tuning"). Another issue in need of investigation would be the dependence of flow splitting on IB canal-position along the rack. According to the recent flume study [9], the splitting coefficient (C) would increase by as much as 25% at locations closer to the inlet.
Although much could be learned from the CFD simulations of flows moving through the rectangular cuboid cavity of Fig 3, the logical next step would be the use of a more realisticallyshaped oral cavity, i.e., curved both laterally and dorsoventrally (Fig 1), and incorporating the varying length of baleen along a rack. This is a challenging task that would require not only an accurate characterization of the oral cavity shape in all three dimensions, but also of baleen chord-wise curvature which changes along the span [38]. Curvature may be important not only for understanding the drag production by those hydrofoils, but also the lift that they may also produce (more energy dissipation) and the possible baleen-tip bending and movement that could follow. The latter would help further understand the wear patterns of baleen [39]. Finally, more work will be needed to clarify the role of the tongue's shape and bulk within the oral apparatus. Per Fig 2, the volume of the tongue is important for the determination of D in and ultimately for drag control.
Modeling Details 1. Input data entry sequence BHC simulations manage information in an order that is opposite to CFD simulations. With the latter, the flows and pressures are calculated first via the numerical solution of the Navier- Stokes equation [26], using input about the morphology and geometry of the body and computational domain. The computation of the drag force follows in a last step, based on the integration of the calculated pressures and velocity gradients along the surface of the object [26]. In the BHC, morphology is also entered first but is followed by the calculation of the pressure P ext external to the PO using an empirically-suggested flow speed. This enables the calculation of the drag force via Eq 3. Calculation of the internal flow speeds becomes the third step of this process, using the solutions of Tables 5-7. The final step involves adjusting and calculating the values of the friction coefficients which then yield the internal pressures. This last step aims at insuring that the BHC-calculated pressures and dissipated viscous energy in the baleen system are consistent with the prescribed flow speeds and predictions of Eq 3. These essential elements are described here and in the five sections that follow. This section begins with input entry, followed by Modeling Details 2 which describes the establishment of the constraints applied in the simulations. The calculations of the flows, pressures and drag are then explained next. The last section, Modeling Details 7, deals with the generalization of the equations listed in Tables 5-7 -valid for a four-baleen-per-rack case-to one with an arbitrary number of baleen (N b ).
The BHC model is simple enough to run in Excel spreadsheet format on any laptop or desktop platform. A simulation begins with the entry of the following basic parameters: the known relevant morphology, namely baleen mean chord (c baleen ), (mean) half-height (h HT ) and separation (d baleen ) (Table 1); the (effective) lateral separation (D in ) of baleen from the tongue wall (Fig 2); fluid characteristics such as seawater mass density and kinematic viscosity, along with the sea level atmospheric pressure (Table 2); and whale swimming speed (U whale ), as well as that of the flows entering the mouth (U in ). (These have been approximated as being the same herein). The number of baleen per rack (N b ) is hard-coded in the spreadsheet as a determinant of the number of cells per column. Working with different values of N b is possible but requires spreadsheet re-configuration.
The value of the empirical flow-splitting coefficient (C = 0.30, 0.44 or 0.53 [9]) is calculated next, via the entry of input a defined as the ratio U Nb /U in of the APT flow speed near the posterior-end of the rack (next to oropharynx), over that of the mouth's entrance. The relationship between input a and C is given by Eq 21b (below).
Next in the sequence is the input of the mean value of the Posterior Opening width (w PO ), calculated from the constraint w PO = D in /1.15 (Table 4) originating from the specifics of the flows external to the body and moving near and outside the PO. These external flows are at the source of the low pressure P ext necessary to maintain the internal flows along and across baleen [3,4]. How the value of P ext and factor "1.15" are arrived at is further explained in Modeling Details 2.
Two final steps are carried out to finalize the friction coefficient profiles (k canal ) in the IB and APT canals (Eqs 13-16). The first consists in adjusting the values of parameter f APT in Eq 15 to yield, via Eq 10, the internal pressure at the PO (P PO ). The latter is usually set at 10-20 Pa higher than the external pressure P ext to insure a favorable pressure gradient just outside the PO. Representative values are shown in Fig 15. (The variations originate from the acceptable values of f APT that put the difference P PO −P ext within this desired 10-to-20Pa interval). The second step adjusts the values of parameters x IB , y rostr , z rostr (Eq 13), to insure a monotonically decreasing pressure on the labial side of each baleen canal (P pi ), i.e., from rostrum to oropharynx (Fig 4). These adjustments are also made so to yield a value of the energy dissipated through all IB canals, and corresponding mouth drag contribution, that is identical to that of the second term in Eq 3. Example values are also shown in Fig 15. Ultimately, both sets of friction coefficient adjustments are necessary to insure consistency with the BHC flow speed solutions (Table 6) which exclude the possibility of backflows into either IB or AP canals. In both frames the circles correspond to a 10m body size (D in = 0.50m (black), 0.30m (red) and 0.72m (yellow)); squares to 15m-case A and triangles for 15m-case B; and starburst for 8m. In all of these cases values of y rostr and z rostr have been set to y rostr = 0.60 to 0.70 and z rostr = 0.0085.

Modeling Details 2. Constraints on the flows external to the posterior opening
A filter-feeding balaenid whale deflects a significant mass of seawater around its head even with the mouth open (Fig 1), effectively accelerating this flow in the vicinity of the PO (U ext ). Such an effect is also seen with the flows moving over the upper surface of any airfoils, typically yielding velocities that are 15%-25% higher than in the freestream [40,41]. These speeds become associated with a zone of low external pressure (P ext ) near the PO, which in turns generates the internal pressure gradients of CFF [3,4]. As such, these external currents interact with the internal flows coming out of the PO (U PO ). The values of U PO , U ext , P ext and P PO are calculated during the input entry process as follows.
Values of U PO versus U ext . The external flow speed is set to U ext = 1.25 U whale per airfoil aerodynamics [40,41]. (Using different values of this coefficient has been considered as well, as further discussed below). On the other hand, the internal flow emerging from the PO is also expected to flow at commensurate speeds, that is, U PO~Uext . On the other hand, the viscous friction within the oral apparatus is expected to reduce the value of U PO somewhat, in comparison to what could be obtained in a world free of viscosity [37]. In nozzles, such an effect is measured through the so-called discharge coefficient, which in the case of open-ended tubes feature outlet velocities varying in the range of 0.7-0.9 of the theoretical, zero-viscosity velocity [37]. In the case of the gently narrowing tubes that approximate the shape of inflating apexvented parachutes (early stage), the associated discharge coefficient stand at about 90 to 94% per a recent CFD study [42]. Approximating each half of the oral cavity as one such tube would yield U PO~0 . 92 U ext~0 .92 x 1.25 U whale = 1.15 U whale . Most importantly, fixing U PO and P PO (P ext ) this way determines the essential physics, that is, the level of viscous friction that in turns yields the drag calculated by Eq 3. With U PO thus assumed, the conservation of the flow mass rate (w PO h HT U PO = D in h HT U in ) is used to arrive at a constraint on the PO width, namely, w PO = D in /1.15 ( Table 4).
Values of P PO versus P ext . Parameter P ext is calculated during the data entry process via Bernoulli's principle, The pressure P in is the known hydrostatic pressure at the mouth and at depth (h HT ) (also calculated during input entry), U in the inlet flow speed (= U whale ), and U ext is constrained to U ext = 1.25 U whale as previously discussed. On the other hand, insuring P PO as greater P ext by about 10-20 Pa is done through the adjustment of the APT friction coefficient profile (see also Modeling Details 4). Knowing the value of U PO , U ext , P ext and P PO basically fix the output of Eq 3's first term. The second term of Eq 3 is also calculated during input data entry, and is used to adjust the values of the IB canal friction coefficient and internal pressures.

Modeling Details 3. Baleen Hydraulic Circuit (BHC) model (four baleen plates)
Flows within oral apparatus are simulated with a low-spatial resolution scheme, illustrated here with a pedagogical 4-baleen example (Fig 3 and Tables 4-6). The more general case of N b baleen plates is discussed further below (in Modeling Details 7) and follows straightforwardly from this example.
The BHC model calculates pressure and flow speeds averaged across pipe/canal sections at the 20 "stations" shown in Fig 3. Because the APT flow loses mass to the baleen canals adjacent to it, a first set of equations is obtained from mass rate conservation, for example as applied to the five-outlet pipe system discussed in [9]: The specific form used along a four-baleen rack is listed in Table 5. Another set of equations is required to characterize how much of the APT flow at a given station (e.g., station i) enters the next IB channel (at station (i+1)): See also Table 5. Generally the flow-splitting coefficient C depends on the location of specific IB canals along the rack as a result of the pressure gradients near the canal's inlet and outlet (see Eq 14 below). But with the intra-baleen canal spacing being small relative to rack length, coefficient C is hypothesized as being only weakly dependent on location, and meaningfully obtained from flow tank data [9]. Note that the station next to the oropharynx (i.e., i = 4 and pi = p5 in Fig 3) involve different flow equations given their location at the end of the baleen rack.
Solving the BHC equations begins by expressing the IB flow speeds in terms of APT speeds via the merging of Eqs 8 and 9. The resulting equations are used then, station by station, to obtain the APT flow speeds, starting with U 1 , then U 2 , and so on. IB flow speeds (U p1 , U p2 , etc.) are then calculated from the known APT speeds (again via the merged Eqs 8 and 9), to yield Tables 5 and 6. The calculation of the pressure drops then follows via Eq 9 (below) using the now known values for U i , U pi (and k canal ; see Eqs 13 and 14 below), resulting in in the equations listed in Table 7.
It is also straightforward to calculate the flow speeds in the APL canal (U Li ) once U i , and U pi have been computed via the flow mass rate conservation (w PO h HT U Ln = d baleen h HT U Pn + w PO h HT U Ln+1 ). An example result is shown in Fig 4.

Modeling Details 4. Pressure drop equations and boundary layer thickness
The pressure profile along the axial length of constant-diameter channels, canals or pipes depends mainly on the size of the boundary layer lining their walls [26,37] and can be calculated from the following relationship relating the pressure drop at the canal's ends to flow speed at its inlet [26]: Typically, U canal is the average flow speed in cases where the pipe flow accelerates, or as the flow speed at the inlet where it does not accelerate [26,37]. The thicker the boundary layer and ensuing energy loss of the flow, the higher k canal and requisite pressure gradient. Because of the assumed flat wall geometry of each canal, the variation of k canal along the baleen rack can actually be partly modelled from the theory of boundary layers over flat plates as described below [26]. In a BHC simulation the pressure drops are calculated from Eq 10, with the flow speeds and friction coefficients calculated in previous steps. How values of k canal are obtained is explained next. Derivation of the coefficient k canal applying to the IB canals ( k IB ) is based on a model of undeveloped pipe flows in which the boundary layer thickness is significantly smaller than the pipe diameter, a result of the larger Reynolds numbers at play ( Table 3). The canal flow moving "outside" the boundary layer is seen as laminar and describable by Bernoulli's principle [26].
(Such a view of Bernoulli flow moving "on top" of the boundary layer is enabled by the concept of the displacement layer [26]). Within the boundary layer, on the other hand, the flow is regarded as turbulent (again, because of the Reynolds number) and following the empirical flow velocity profile of along a flat plate [26], which grows over the distance x as follows (measured from the inlet): With the kinematic viscosity ν ranging from 1.05 to 1.45 x10 -6 m 2 /s, one typically finds small boundary layer thicknesses at the outlet of each IB canals, namely, δ~0.003m with the baleen channels located near the rostrum (where U in = 0.8 m/s), and δ~0.0048m near the oropharynx (where U in~0 .06 m/s). These values are indeed small in comparison with the 0.01mbaleen baleen separation defining the width of IB canal (which are typically 0.07m-long laterally). With undeveloped flows, the presence of the (displacement) boundary layer growing downstream creates an ever-smaller constriction for the Bernoulli flow to pass through, thereby accelerating it and reducing its ambient pressure [26]. Such gain in speed at the outlet is expressed in terms of the inlet speed via conservation of the mass flow rate through the canal as follows [26]: The length δ outlet corresponds to the boundary layer thickness at the IB canal outlet, as evaluated from Eq 11 (with x = c baleen ). The factor 1/8 arises from the calculation being performed with respect to the displacement thickness, a length scale that is smaller than δ (Eq 11) by a factor 1/8 [26]. The pressure drop along both ends of a given IB canal ( ΔP i ) can be determined by merging Eqs 10 and 12 with U inlet = U pi and Bernoulli's equation , to yield k IB (i) = (d baleen /(d baleen −δ(i)/4)) 2 -1 [26]. Here the index i represents the baleen canal being considered along each rack. (The distance to the rostrum is given by i d baleen ).
To this value of k IB one must add (first) the friction losses attributed to the two bending paths a given IB current has to make to ultimately reach the PO, namely, from APT to IB, and then from IB to APL. In hydraulics, these are attributed to so-called secondary flows which move fluid off-axis within each bend [26,35,36]. Given the associated values of the Reynolds numbers (~10 3 ; Table 3), such losses significantly rises the value of k IB , namely up to 1.0 to 3.0 per units of ΔP/0.5ρU 2 [35,36]. The IB canal friction coefficient that results is written as follows: k canal k IB ðiÞ ¼ x IBbend f1 À y rostr e À z rostr i g þ d baleen The contribution of the two bends is represented by coefficient x IB which is set in the range of 0.5 to 2.5 (Fig 15), following a comparison of the second term in Eq 3 with a direct calculation of energy dissipation and drag produced (Eq 20) . Such values turn out as quite typical in the microfluidics of bends [35] and hydraulics of small diameter pipes [36]. The function appearing within the curly brackets in Eq 13 is an ad hoc addition designed to enable a monotonically-decreasing pressure along the APL canal (P pi ), typically in the anterior third portion of each rack (Fig 4). The value of y rostr and z rostr deemed necessary to achieve this are in the range y rostr = 0.60-0.70 and z rostr = 0.0085. . The latter indeed turns out smaller than half the canal width as expected (undeveloped flows) and over most of a rack. Where the boundary layers meet, the corresponding k IB -values are larger-but by only a factor 1.05x -than the values found in pipe systems supporting fully-developed and turbulent flows at the same Reynolds number (~10 3 -10 4 ) [26].
Modeling the friction coefficient in the APT canal proceeds similarly. With this canal formed by the tongue and baleen rack itself (Fig 1), there is only one continuous wall which in turns supports only one continuous boundary layer. Further, loss of fluid mass through the adjacent IB canals results in decreasing APT flow speeds per mass rate conservation [9] (Fig  4). Modeling APT pressure drops at each station again begins with using Eq 10 but now applied to "pipe" sections connected in series, each of length spanning baleen separation d baleen . These pipes are again assumed to support undeveloped flows, with a Bernoulli flow accelerating atop of the displacement boundary layer growing along the lingual wall. Per the mass rate conservation one has U inlet h HT D in = U outlet h HT (D in -δ (outlet) /8) (with U inlet = U i ), which in turns yields: The effect of adjacent IB canal on overall friction is accounted through f APT , which is tuned during input-entry so to lower P Nb and P pNb to values slightly above P ext (Fig 4; and Modeling Details 2), resulting in the values shown in Fig 15. An example of boundary layer thickness along the APT canal calculated via Eq 16 is also shown in Fig 16. The algorithm based on Eq 10 yields pressure profiles along the APT and APL canal which are continuous within numerical discretization variability (Fig 4). It follows that they yield the same pressure difference between any two points in the cavity, via integration of pressureincrements along arbitrary paths in the APT and/or APL canals starting and ending at those same two points. For example, calculating the pressure difference (P 6 -P p6 ) directly from Eqs 30 and 32 yields a value of 51.94Pa, in comparison to 51.78Pa obtained from integrating the pressure decrement (P i+1 -P i ) in the APT canal from site i = 6 and to site i = 200, plus the difference P 200 -P p200 (following a path through IB canal #200), and finally adding the integration of the decrements (-P pi+1 + P pi ) in the APL canal from i = 200 to i = 6. The 0.3% difference that results is typical of the many paths that have been tested.

Modeling Details 5. Viscous friction drag of the oral apparatus (Eq 3)
The drag contributed by the oral apparatus is obtained by considering the rate of energy dissipated by the flows moving through it. This follows from the fact that with bodies traveling through a viscous fluid at constant speed and under the action of a single external and forward-directed force (thrust), the work done by the latter turns out as identical to the energy (E friction ) dissipated by the viscous friction sustained by, and turbulence imparted to the fluid. With drag being equal to thrust, and over travel time scale T, one therefore has F D = E friction /(U whale T). (Note that this approach would not be correct if the body were to accelerate, since drag is no longer equal to thrust, and more importantly, since drag generation would then consume more energy beyond viscous friction and turbulence, i.e., in the form of the extra kinetic energy gained by the accelerating fluid. This latter effect leads to the socalled acceleration reaction or added mass drag [10]).
In what follows, scale T will be chosen as the time needed for the whale to swim a distance equal to one body length. Eq 3 is derived by estimating the rate of energy dissipation (E friction /T) inside the mouth, here seen as a juxtaposed pair of boxes each featuring laterally-offset inlets and outlets, and carrying a rack of hydrofoils (baleen) of known drag characteristics. In this model, the energy dissipated anteroposteriorly (AP) is calculated independently from that dissipated by the lateral (IB) flows moving past the hydrofoils.
Generally, the rate of viscous dissipation arises from a turbulence component, dominant at high Reynolds numbers, superposed to a frictional component that dominates at low Reynolds number. The effects of the former can be isolated following an argument based on Reynoldsaveraging [43], wherein the flow field is decomposed into a time-averaged flow speed U and a time-dependent fluctuating component u for which the time-averaged values are u ¼ 0 and u 2 6 ¼ 0. With a 2-dimensional velocity field (i.e., longitudinal and lateral) going through any  through the baleen/hydrofoil rack matching the mass flow rate thought the mouth's inlet: Parameter L mouth corresponds to the length of the oral apparatus, and factor ½ (D in + w PO ) serves as a travel distance scale for the fluid moving laterally within the cavity. The latter is adjusted by the quotient appearing in between square brackets to ensure conservation of the mass flow rate everywhere. Within the assumptions of the BHC, this length scale and V lat basically grow with mouth inlet size as D in /L mouth (in accord with the solution shown in Table 6). The rate of dissipated energy for the entire cavity is thus obtained by adding to the energy dissipated longitudinally (in the right-hand-side of Eq 17b but without the U lat -terms) the energy lost thought the lateral motions through the hydrofoil racks. Multiplying those energy densities by the volumetric rate D in h HT U in yields the following drag force generated by open mouth: The factor "2" multiplying the large curly brackets adds the effects due to both baleen racks; and setting U in as equal to U whale leads to Eq 3.
Calculating the dissipated energy incurred by the lateral flows though baleen/hydrofoils also yields an expression of the body drag contributed by the dissipated energy in each of the 2N b canals (= k IB ½ ρ w U pi 2 ), namely This result follows from using Eq 17b in each IB canal, associating "in" with the inlet and "PO" with the outlet, setting the "lat" speeds to zero, and having the same "long" speeds at both inlet and outlet. From this one obtains e friction long = P in −P PO , which along with Eqs 10, 13 and 14 yields the energy density (J/m 3 ). Multiplying the latter by volumetric rate d baleen h HT U pi gives the rate of energy dissipation (Watts) in a canal. Summing over all baleen canals and multiplying by the time T spent by the whale to travel a distance equal to its body length (L body ) yields the work (F D mouth | IB L body ) needed from the whale to replace the energy dissipated in the IB canals and ensuing mouth drag. Eq 20 is another important element of the BHC as it allows, via the comparison with the second term of Eqs 3 or 19, the determination of the final value of k IB via the tuning of x IB .
Modeling Details 6. APT flow speeds near the oropharynx and relation to body size The AP speed at the end of the APT canal is interesting as it corresponds to the speed of the ready-to-eat slurry near the oropharyngeal wall. One can calculate it's ratio to mouth inlet flow speed as follows (in a 4-baleen system; see Table 6): Cross section areas A IBcanal and A in are those of the IB and APT canals respectively (Fig 2 and  Table 4) ( Table 4). In the case of a more general N b = N-1 baleen system this ratio would be calculated as (see Modeling Details 7): In the limit where the number of baleen plates (and N) becomes very large, the value of a converges to a near-zero value since the ratio C A IBcanal /A in is less than unity. Note also that this function is also very close numerically to the exponential a~exp(-N b A baleen /A in ) in the limit of large baleen plate numbers.

Modeling Details 7. BHC modeling for N b baleen (per rack)
Generalizing the BHC approach to an arbitrary number N b of baleen per rack (= N-1) is straightforward once the 4-baleen case is understood. One begins with the mass rate equation which insures that, over a given time interval, the fluid mass entering a selected APT canal ( Fig  3) is equal to the fluid mass exiting through both of the canal's IB channel and APT outlet: Again both cross section areas A IBchannel and A in known (Table 4). Note the equation for the last pipe section where there is no APT channel to discharge into, per Fig 3 (assuming a closed oropharynx). Eq 22 ensure that the mass rate entering the tongue-baleen gap canal (A in U in ) is identical to the flow rate leaving each posterior opening (PO): This can be shown by noting first that the flow speeds (U Li ) in the baleen-lip canal (APL canal) arise from the merging of the flows exiting each IB canal, namely, and from the flow mass rate through the PO is given by Note that the flows U Li defined above exclude the flow U L0 that enters the APL canal at the rostrum, hence the extra term on both sides of the mass conservation equation above.
The next series of equations tracks the flow splitting: As discussed in the text, the value of the flow splitting coefficient C is obtained from flow tank data [9] and is assumed to be the same along the full length of each baleen rack. Merging the above equations and solving for flow speeds proceeds hence, beginning at station n = 1, then solving at station 2, and so on. Along the sections of the APT canal the solution is as follows: On the other hand, each IB channel features a solution of the form: The flow speeds in the APL channel are obtained by using Eqs 27 and 28 along with mass conservation enforced in this canal (Fig 3): With the flow speeds known everywhere in the buccal cavity, one then calculates the values of the friction coefficients k IB and k APT from Eqs 13-16. The corresponding pressure drops then follow via Eq 10: P nÀ 1 À P n ¼ k APT ðnÞ 1 2 r w U n 2 1 n ðN À 1Þ ð30Þ P n À P pn ¼ k IB ðnÞ 1 2 r w U pn 2 1 n ðN À 1Þ ð31Þ Once the flow speeds and pressure drops in both APT and IB canals are computed, the friction coefficients k APL in the APL canal is reconstructed using the associated Darcy-Weisbach equation (Eq 10): k APL ðnÞ ¼ P pnÀ 1 À P pn 1 = 2 r w U Ln