## ABSTRACT

Countless aquatic animals rotate appendages through the water, yet fluid forces are typically modeled with translational motion. To elucidate the hydrodynamics of rotation, we analyzed the raptorial appendages of mantis shrimp (Stomatopoda) using a combination of flume experiments, mathematical modeling and phylogenetic comparative analyses. We found that computationally efficient blade-element models offered an accurate first-order approximation of drag, when compared with a more elaborate computational fluid-dynamic model. Taking advantage of this efficiency, we compared the hydrodynamics of the raptorial appendage in different species, including a newly measured spearing species, *Coronis scolopendra*. The ultrafast appendages of a smasher species (*Odontodactylus scyllarus*) were an order of magnitude smaller, yet experienced values of drag-induced torque similar to those of a spearing species (*Lysiosquillina maculata*). The dactyl, a stabbing segment that can be opened at the distal end of the appendage, generated substantial additional drag in the smasher, but not in the spearer, which uses the segment to capture evasive prey. Phylogenetic comparative analyses revealed that larger mantis shrimp species strike more slowly, regardless of whether they smash or spear their prey. In summary, drag was minimally affected by shape, whereas size, speed and dactyl orientation dominated and differentiated the hydrodynamic forces across species and sizes. This study demonstrates the utility of simple mathematical modeling for comparative analyses and illustrates the multi-faceted consequences of drag during the evolutionary diversification of rotating appendages.

## INTRODUCTION

The diverse array of aquatic animals that rotate their appendages for locomotion and prey capture offers rich material for considering the role of hydrodynamics in morphological and kinematic diversification (Fish, 1984; Johansson and Lauder, 2004; Koehl, 1996; McHenry et al., 2003; Ngo and McHenry, 2014; Richards, 2010; Webb and Blake, 1985). Although the drag forces in terrestrial systems are often negligible, the aquatic appendages of dragonfly larvae, snapping shrimp and mantis shrimp, for example, necessarily incur substantial drag (Anker et al., 2006; McHenry et al., 2012; Tanaka and Hisada, 1980; Versluis et al., 2000). Therefore, relationships among the parameters that influence drag have the potential to inform interpretations of morphological diversity.

*F*

_{d}) is commonly modeled with the following equation:

where ρ is the density of the fluid medium, *S* is a characteristic surface area of the object moving through the fluid, *v* is the linear velocity of the object relative to the fluid and *C*_{d} is the drag coefficient, a non-dimensional measure of shape. The drag coefficient is typically determined from drag measured by a force transducer for a body exposed to the uniform flow generated by a flume (e.g. Van Wassenbergh et al., 2015b). In contrast, a rotating body is exposed to flow that varies linearly along its length. A blade-element model of drag may account for this position dependency as the sum of elements along the appendage's length (Blake, 1979). For each instant in the rotation of an appendage, this sum may be formulated from measurements of the dimensions of each element, probable values for its drag coefficient, and its position-dependent velocity (McHenry et al., 2003). This approach treats all elements as independent force-generating units and thereby neglects factors such as span-wise flow and the high shear stress that may be generated at the distal end (e.g. Dickinson et al., 1999; Lentink et al., 2009). The relative simplicity of the blade-element approach yields a method for calculating drag that is computationally efficient and therefore amenable to comparative analysis (Richards, 2010; Walker, 2004). However, it is unclear whether the simplifying assumptions of a blade-element approach are violated for any particular rotating appendage. Therefore, the present study compared the predictions by blade-element modeling against both experimental approaches and more elaborate computational analyses.

To fully incorporate the effects of rotation on hydrodynamics, many studies have employed computational fluid dynamics (CFD) (e.g. Jiang and Kiørboe, 2011; Li et al., 2012; Liu et al., 1996; Van Wassenbergh et al., 2015a). CFD analyses can consider the fluid forces generated by complex geometries with fully resolved volumetric flow fields. However, the investment in the development and computational resources to perform CFD simulations are sufficiently intensive that they typically preclude large-scale comparative, evolutionary analyses that incorporate detailed morphological variation. An analysis of rapidly rotating pipefish snouts found that CFD models agreed with the predictions of a more simple blade-element model, where drag was treated as the sum of independent structural elements with linear flow assumed (Van Wassenbergh and Aerts, 2008). Similarly, an analysis of a hovering fruit fly also yielded convergent results between blade-element and CFD models (Walker, 2002). In contrast, in the context of the more complex aerodynamics of flapping flight, a modified CFD analysis performed better than a blade-element approach (Nakata et al., 2015).

*c*chord length; linear dimension in direction of flow (m)

*C*_{d}drag coefficient (dimensionless)

- CFD
computational fluid dynamics

*E*_{d}energy lost due to drag (J)

*F*_{d}drag (N)

*h*width; linear dimension perpendicular to flow (m)

*i*element number

*k*shape coefficient (dimensionless)

*l*thickness; linear dimension along the longitudinal axis of the structure (m)

*L*characteristic length (m)

- micro-CT
micro-computed tomography

- MSE
mean-squared error

*n*number of elements

- PGLS
phylogenetic generalized least squares

*r*distance to pivot point (m)

*Re*Reynolds number (dimensionless)

*S*surface area (m

^{2})*t*time (s)

*T*strike duration (s)

*T*_{d}drag-torque index (dimensionless)

*U*velocity of fluid in flume (m s

^{−1})*v*velocity of object relative to fluid (m s

^{−1})- γ
angle of dactyl/propodus rotation (rad)

- μ
dynamic viscosity (Pa s)

- ρ
density (kg m

^{−3}) - τ
drag torque (N)

Mantis shrimp (Stomatopoda) offer an intriguing system in which to examine the effects of shape, size and velocity on drag forces (Fig. 1). ‘Smasher’ mantis shrimp evolved from ‘spearers’ (Porter et al., 2010) and, during this transition, the dactyl-open evasive prey capture motion characteristic of spearers switched to a primarily closed-dactyl, hammering motion directed at hard-shelled prey (Patek, 2015). Along with this shift in the orientation and target of the raptorial strike, smashers dramatically increased in speeds and accelerations (Table 1) (Cox et al., 2014; deVries et al., 2012; Kagaya and Patek, 2016), evolved enhanced elastic energy storage capabilities (Patek et al., 2013; Rosario and Patek, 2015), yet also experienced a substantial decrease in range of body size and appendage size compared with their ancestral spearers (Anderson et al., 2014; Blanco and Patek, 2014). Thus, the evolutionary shift to smashing encapsulates the major axes of drag: shape, speed and size.

The goal of the present study was to analyze the hydrodynamic consequences of kinematic and morphological variation in the fast rotational motion of the mantis shrimp raptorial appendage. We evaluated the methods of drag measurement, blade-element modeling and CFD to quantify drag on this rapidly rotating structure. In addition, we performed a new kinematic analysis of a spearing mantis shrimp, *Coronis scolopendra*, and conducted phylogenetic comparative analyses of the relationships among drag-related parameters across mantis shrimp. We considered the torque generated by drag (henceforth referred to as ‘drag-torque’) and the propensity of an appendage to generate this torque, which was measured with the drag-torque index (McHenry et al., 2012). These approaches allowed us to address three guiding questions: (1) how does the balance of evolutionary variation in kinematics, shape and size influence variation in drag across stomatopods; (2) given the wide range of high Reynolds numbers (*Re*) for this system (Table 1), can we make simplifying assumptions about drag calculations, in the context of the added complexities of rotational, rather than translational, motion; and (3) given the results of the comparative analyses and fluid dynamic methodological comparisons, which aspects of mantis shrimp fluid dynamics are most relevant to evolutionary diversification and should be the primary focus of future fluid dynamic analyses in mantis shrimp?

## MATERIALS AND METHODS

### Strike kinematics of *Coronis scolopendra*

We supplemented previously published kinematics of very small and larger spearers through high-speed videos of a medium-sized spearer, *Coronis scolopendra* Latreille 1828 (Crustacea: Stomatopoda: Lysiosquilloidea: Nannosquillidae). Individuals were collected in Florida (Florida Division of Marine Fisheries Management, license no. SAL-13-1278-SRP) and kept in aquaria (24–28°C, 34–36 ppt salinity), where the animals formed their own burrows (sugar-sized Oolite, Aragonite, CaribSea, White City, FL, USA). They struck live brine shrimp introduced by pipette or forceps and were filmed using a digital high-speed imaging system (10,000–15,000 frames s^{−1}, 1/10,000–1/15,000 shutter, 512×512–512×256 pixels, APX-RS, Photron Inc., San Diego, CA, USA). A ruler was placed in the plane of the striking animal to calibrate each series of videos.

Five points were tracked manually in each image sequence of a raptorial appendage strike (Fig. 2). Tracking began when the carpus rotation started, and ended when the appendage made contact with the prey or feeding device (MATLAB v. R2011a, The MathWorks, Natick, MA, USA). The five points included a central point on the merus as well as points on the distal ends of the merus, carpus, propodus and dactyl (Fig. 2). The percent error in each measured displacement averaged 9.2±3.5% (range 3.5–16.5%), assessed by 10 repeated measures of a representative strike.

Kinematics were calculated from the coordinate measurements (MATLAB v. R2013b, The MathWorks). Rates of change were calculated from the first and second derivatives of a least-squares fit to a 10th-order polynomial determined to positional data (Cox et al., 2014; deVries et al., 2012). If maximum speed occurred at the beginning or end of the digitized portion of the strike, then the strike was eliminated from the dataset, because the true peak speed could have occurred outside the digitized segment of the strike. Using linear models, we tested whether kinematics were correlated with body size (R v. 3.0.2, R Foundation for Statistical Computing, Vienna, Austria) using propodus length as a proxy for body size (Claverie et al., 2011; deVries et al., 2012; Patek and Caldwell, 2005).

### Strike kinematics across species

We used previously published and *C. scolopendra* strike kinematics to determine scaling relationships. We incorporated kinematics from *Lysiosquillina maculata* and *Alachosquilla vicina* (deVries et al., 2012), *Gonodactylus smithii* (Cox et al., 2014; Patek et al., 2007) and *Neogonodactylus bredini* (Kagaya and Patek, 2016) in a form consistent with present measurements. We measured the propodus length (length of the propodus between the propodus–dactyl joint to the proximal-most point visible on the propodus) and the striking body length (length from the propodus–dactyl joint to the insertion point of the lateral extensor muscle on the dorsal surface of the carpus) from archived digital images (Fig. 3).

We performed phylogenetic generalized least squares (PGLS) regression on logged species averages against logged average striking body length to examine the relationship between size and kinematics across species. Maximum kinematic values for each individual were used to calculate the average maximum value for each species. We incorporated a pruned and time-calibrated molecular phylogenetic tree based on nucleotide sequence data (Porter et al., 2010) with hard-bound calibration points from fossil data to establish time calibrations (Claverie and Patek, 2013). PGLS analyses were performed in the R package ‘caper’ (https://cran.r-project.org/web/packages/caper/index.html), with delta and kappa fixed at unity and lambda estimated using maximum-likelihood methods, which allows the model to deviate from a strict Brownian motion.

*Gonodactylus smithii*was employed as a representative smasher,

*L. maculata*represented large spearers and

*A. vicina*served as a characteristic small spearer (Fig. 3). We applied a sensitivity analysis that examined the effects of kinematics on the maximum torque and energetic expense generated by drag (described below). The striking motion, defined as the angle of the propodus–dactyl unit during a strike with respect to time, approximated a sigmoidal curve for all species. The species primarily differed in terms of strike duration and angular excursion of the appendage. We found that the excursion, γ

_{d}(in radians), was correlated with the strike duration (in seconds),

*T*, as specified by the following equation:

Therefore, much of the variation in strike kinematics across species may be explained by differences in strike duration.

*L. maculata*that exhibited the characteristic sigmoidal pattern that we found in all species (Fig. S1). After normalizing for the strike duration and maximum excursion, we used a non-linear least-squares fit to characterize this strike with the following fifth-order polynomial:

where and are the respective values of normalized angular position and time. This pattern was used to vary the kinematics across species by multiplying normalized time values by the strike duration and the normalized angular position by the corresponding angular excursion (using Eqn 2). This approximation was validated by comparing the performance for kinematics of Eqn 3 against the measured kinematics of equal duration using a Kolmogorov–Smirnov goodness-of-fit test for *G. smithii* (drag energy: *P*=0.64; maximum torque: *P*=0.47; *N*=47), *C. scolopendra* (drag energy: *P*=0.36; maximum torque: *P*=0.68; *N=*50) and *L. maculata* (drag energy: *P*=0.16; maximum torque: *P*=0.08; *N=*58). Calculations of the mean-squared error (MSE) for the fit (*G. smithii*: MSE=0.039 rad^{2}; *C. scolopendra*: MSE=0.005 rad^{2}; *L. maculata*: MSE=0.262 rad^{2}) indicated that the spearer, *L. maculata*, exhibited substantial variation not represented by Eqn 3.

### Drag measurements on physical models

where *U* is the velocity of flow, *L* is a characteristic length, and ρ and µ are, respectively, the density and dynamic viscosity of water. The characteristic length is commonly selected as a dimension in the direction of flow, whereas we used the propodus length because of the irregularity of the cross-sectional shape in the direction of flow. The 3D prints were enlarged to compensate for the slow flow speeds generated by the flume. Prints for two of the taxa did not have published kinematic data, so the scaling was based on speed from similar taxa. In particular, *Gonodactylaceus* *falcatus* models were assumed to have speed similar to *G. smithii*, and *Hemisquilla* *californiensis* were approximated by the speed of *L. maculata*. For the three relatively slow species (*H. californiensis*: undifferentiated; *L. maculata* and *C. scolopendra*: spearers), we were able to match and exceed the *Re* of measured strikes (Table S1). We used the smaller end of the size range for the smashers *G. smithii* and *G. falcatus* to correctly match their *Re* while accommodating the limits to the size of the 3D printer (Z-Corp 310, 3D Systems, Rock Hill, SC, USA) and the top speed of the flume (1 m s^{−1}, model 2436, Rolling Hills Research, El Segundo, CA, USA) (Table S1).

The 3D printed models were based on scans of the dactyl and propodus segments acquired from micro-computed tomography (micro-CT; Model HMXST225, X-Tek, Nikon Metrology NV, Leuven, Belgium). The carpus, which is proximal to the propodus, was excluded because of its modest contribution to the generation of drag (McHenry et al., 2012). The scans were performed for two smasher species (*G. smithii* and *G. falcatus*), two spearer species (*C. scolopendra* and *L. maculata*) and an undifferentiated species (*H. californiensis*) (Figs 1, 3). The individual dactyl and propodus segments were digitally separated (Mimics v13.0 for X64, Materialize, Belgium) and their wetted surfaces were identified and smoothed (Geomagic Studio 11, 64-bit edition, Geomagic, Cary, NC, USA). The resulting geometries were printed and then hardened with cyanoacrylate to create enlarged physical models of the dactyl and propodus that could be mounted with the dactyl in either an opened or closed position relative to the propodus (Fig. 1).

We measured the drag on the physical models over a range of flow. Each physical model was suspended in the working section of the flume with the long-axis of the propodus perpendicular and the dactyl opened toward oncoming flow (Fig. 4). The sting holding each model was positioned such that force was applied in the direction of flow velocity (Fig. 4) against a force sensor (DS2-4, Imada, Northbrook, IL, USA) for a 1-min recording (SW-1 Data Acquisition package, Imada). The drag generated by the sting was subtracted from each measurement.

### Computational fluid dynamics

We tested a blade-element model of the hydrodynamics of a rotating structure with CFD. We animated the CFD model with the strike kinematics of the smasher, *G. smithii*, and approximated its geometry with a series of 19 frustum segments (virtual slices of the appendage) connecting 20 ellipses that match the cross-sectional shape along the length of the appendage (Fig. S2). This multi-frustum model was created in the CAD program GAMBIT v. 2.4.6 (ANSYS, Canonsburg, PA, USA) and then imported into ANSYS DesignModeler 14.5.7, where it was surrounded by a spherical flow domain with a radius of 0.05 m. The flow domain was meshed with 11 million tetrahedrals, with a density decreasing with distance away from the appendage (growth rate 1.1 per layer starting from a mesh element size of 20 µm, ANSYS Meshing 14.5.7). This density was sufficient, as suggested by a simulation with one-third of the number of tetrahedrals, which predicted the maximum hydrodynamic torque to differ by only 1.5%. The model was subjected to a constant acceleration of 4×10^{6} rad s^{−2}, from standstill to 2800 rad s^{−1} in 0.70 ms, which corresponds closely to the values (2821 rad s^{−1} at 0.72 ms) for *G. smithii* (Table 1). Simulations solved the flow field for the full Navier–Stokes equations for unsteady laminar flow (ANSYS Fluent 14.5.7) with a fixed time-step (3.5 µs). Quartering the time-step size had a negligible effect on the calculated torque (<0.7% difference). In these simulations, the mesh of the entire flow domain was rotated at constant acceleration with respect to the fixed reference frame (DEFINE-ZONE-MOTION function in Fluent). The outer boundary of the spherical flow domain was defined as a pressure outlet with zero-gauge pressure and backflow perpendicular to the no-slip boundary of the appendage.

The torque generated by drag was determined by summing the torques due to viscous and pressure forces on each surface element of the multi-frustum model. This calculation was compared with the prediction for the same conditions by a blade-element model, as done previously (Van Wassenbergh et al., 2008). As detailed below, a blade-element model calculates the quasi-steady fluid forces as the sum of forces generated by virtual sections along the structure. In this implementation, the mantis shrimp appendage was divided into 19 segments and the torques generated by drag and added mass were calculated for each interval of time (eqns 2.6 to 2.15 in Van Wassenbergh et al., 2008).

We compared the drag generated by rotational motion with that generated by linear motion. This was achieved by CFD simulations with the multi-frustum model in steady, linear flow of 20 m s^{−1} with a cylindrical flow domain (0.3 m in length, 0.05 m in radius). The appendage was oriented perpendicular to the central axis of the cylinder. Pressure inlets and outlets at the ends of the cylinder induced flow along this axis toward the leading edge of the appendage. The walls of the cylinder were defined with a velocity equal to freestream flow (20 m s^{−1}) and the multi-frustum surface was defined with a no-slip condition. We used Menter's shear stress transport model in ANSYS Fluent (Menter, 1994) to resolve the flow patterns in the wake, which accurately calculates the steady-state drag of bluff bodies at similar *Re* (Goyens et al., 2015; Van Wassenbergh et al., 2015b). A mesh convergence analysis showed that the drag did not vary substantially when refining from 11 million to 16 million cells (+0.1%) (Fig. S3). We retained the latter density for our simulations (mesh size at the multi-frustum surface of 1.6 µm), which showed convergence after 2000 iterations. The drag coefficients of each individual frustum segment, as well as the overall drag on the model, were calculated from the solution.

To evaluate the influence of the drag coefficient (*C*_{d}) on the accuracy of the blade-element model, three versions of blade-element models were compared: (1) a model using *C*_{d} values from measurements of a range of long elliptical cylinders with different aspect ratios (see Eqn 5 below), (2) a model that used the *C*_{d} values calculated for each individual frustum segment of the steady-flow CFD simulation described above, and (3) a model that used the total drag force from the steady-flow CFD simulation. To focus on the effect of drag forces, differences in added-mass torque between the blade-element models and CFD at the instant of pure added mass resistance (the first time steps of the simulation when velocity is very low and drag is negligible) were cancelled out by adding a constant torque (8.54×10^{−5} N m or 19% of the CFD's added mass torque). This analysis allowed us to test the case of an extreme rotational acceleration, as observed in the smashing mantis shrimp.

### Modeling torque from drag measurements

We modeled the torque generated by drag during the strike of a raptorial appendage from our measurements of dynamically scaled models. This necessitated a consideration of the differences in flow between a rotating structure (as in mantis shrimp) and one exposed to the translating flow in a flume (as in our drag measurements). We modeled drag torque using a quasi-steady blade-element approach that treated the morphology as a series of chord-wise (i.e. in the direction of flow) elements that vary in their dimensions to conform to the shape of the dactyl–propodus unit.

*F*

_{d}) was considered equal to the sum of drag on all elements, as indicated by the following equation (Batchelor, 1967):

*l*is the thickness (i.e. linear dimension along the longitudinal axis of the structure) and

*h*is the width (i.e. linear dimension perpendicular to flow) of each element

*i*, with

*n*equaling the total number of elements. This measure of drag coefficient, (

*C*

_{d})

*, may be distinguished from the coefficient for the entire appendage, because it represents the contribution of an individual element. The drag coefficient for each element was modeled in a form similar to a uniform elliptical cylinder, given by the following equation (Hoerner, 1965):*

_{i}where *k* is a shape coefficient that varies with the geometry of the structure and *c* is the chord length (i.e. linear dimension in the direction of flow). We measured the dimensions (*c* and *h*) from our micro-CT scans with the dactyl–propodus unit (with the dactyl closed against the propodus) for 20 evenly spaced elements. For a uniform elliptical cylinder, a fixed value for the shape coefficient (*k*=0.015) is predictive of empirical drag measurements (Hoerner, 1965). The shape coefficient, *k*, was determined for each of our dynamically scaled models from our measurements of drag using nonlinear least-squares (‘lsqcurvefit’ in MATLAB) to minimize deviation from our drag measurements across flow speeds.

*/*d

*t*is the angular velocity and

*L*is length of the dactyl–propodus unit (i.e. equal to the characteristic length for

*Re*, Eqn 4). The drag-torque index,

*T*

_{d}, indicates the propensity of an appendage to generate drag torque during a strike. Its calculation incorporates the variation in

*C*

_{d}(Eqn 6) and associated shape coefficients (

*k*) of blade elements that we determined from drag measurements, as articulated by the following equation (McHenry et al., 2012):

where *r* is the distance between an element and the pivot point of the dactyl–propodus unit, approximated as the proximal-most point on the propodus. We calculated drag torque index, *T*_{d}, across all of our measured species with the dactyl in the open and closed positions.

### Effects of shape, size and kinematics on strike performance

*E*

_{d}) which was calculated by integrating drag torque over the angular displacement of a strike (McHenry et al., 2012):

where γ_{d} is the total angular excursion over a strike. This integral was solved numerically using the ‘trapz’ function in MATLAB with the strike excursion divided into 1000 equal intervals of time.

We used our model of drag torque (Eqn 7) to perform a sensitivity analysis that explored the effects of kinematics, size and shape on maximum drag-torque and drag energy. Each simulation independently modified strike duration, the length of an open or closed dactyl–propodus unit, and drag-torque index over a series of 100 simulated strikes. We varied each parameter of interest by 0.75 orders of magnitude above and below the mean value for the species while setting all other parameter values equal to the mean value for each species. For each simulation, we calculated the maximum torque (Eqn 7) and drag energy (Eqn 9).

## RESULTS

### Strike kinematics

*Coronis scolopendra*'s strikes followed a sequence typical of spearer species (Table 2): (1) the dactyl opened from the propodus, (2) the propodus slid along the merus, (3) the propodus rotated toward the prey and (4) the dactyl and propodus made contact with the prey. The slide of the propodus prior to propodus rotation indicated that the strikes were powered by spring-loaded skeletal elements. Prey distance from the mantis shrimp's eye was on average 7.12 mm (range±1 s.d.=3.68–9.94±1.96 mm). The distal tip of the propodus reached an average peak linear speed of 2.1 m s^{−1} (0.7–3.6±0.6 m s^{−1}) and average peak linear acceleration of 2600 m s^{−2} (300–11,000±1200 m s^{−2}). Peak linear velocity was reached from 0.2 to 6.1 ms after the propodus began rotating, and strike duration ranged from 2.3 to 13.1 ms.

*Coronis scolopendra* size and kinematics were not correlated when analyzed in terms of propodus length and maximum linear velocity (linear regression: *R*^{2}=0.07, *N*=6, *F*_{1,5}=0.377, *P*=0.57). A non-significant trend was present in the analysis of propodus length and maximum linear acceleration (linear regression: *R*^{2}=0.57, *N*=6, *F*_{1,5}=6.74, *P*=0.05).

The phylogenetic (PGLS) regression analyses of *C. scolopendra* and five previously studied species revealed a significant negative association between angular kinematics and striking body length (Fig. 5, Table 1) (angular velocity: *P*=0.02, *R*^{2}=0.76; angular acceleration: *P*=0.02, *R*^{2}=0.77). Linear kinematics were not significantly associated with size, although linear acceleration appeared to follow a non-significant negative relationship with size (*P*=0.06, *R*^{2}=0.6).

### Computational fluid dynamics

We performed CFD simulations to test the accuracy of a blade-element model in estimating the drag on a rotating structure at the scale of a raptorial appendage. The simulation of steady, linear flow of 20 m s^{−1} on the multi-frustum model of the appendage of *G. smithii* showed considerable span-wise flow at the leading edge directed towards the free ends (frustum segments 1 to 3, and 17 to 19; numbered from proximal to distal), coupled with lower calculated drag coefficients at these outer regions compared with literature values for infinitely long elliptical cylinders (Fig. 6A). Span-wise flow at the leading edge in the proximal direction towards the mid-region of the appendage model was observed on the proximal side of the bulging part of the appendage (frustum segments 11 to 15). Drag coefficients for long elliptical cylinder from Eqn 5 for *k*=0.015 consistently overestimated the values calculated by the CFD simulation (frustum segments 10 to 12), but the difference was most prominent for the above-mentioned frustum segments where span-wise flow was observed (Fig. 6A). According to this CFD simulation, total drag force would be overestimated by 55% when using drag coefficients from long elliptical cylinders in a blade-element model.

The simulation of accelerated rotation of the multi-frustum model revealed a zone of high positive pressure that functioned to resist rotation at the leading edge of the elliptical cylinder. This was accompanied by a zone of high negative pressure on the proximal and lateral surfaces of the multi-frustum model (Fig. 6B,C). The performance of a blade-element model (Van Wassenbergh et al., 2008) in approximating the torque about the fixed center of rotation, as calculated by CFD, depended on the *C*_{d} input treatment (Fig. 6D–F). The model using the *C*_{d} value of a long elliptical cylinder (Eqn 5) (Hoerner, 1965) showed the largest difference: a steeper increase in the resistive torque with increasing angular velocity owing to drag forces caused this blade-element model to overestimate the final torque (at time 0.700 ms) by 54% (Fig. 6D). The difference between the two other models and the CFD solution was much smaller: the effect of drag was underestimated at the final simulation time by 4.6% in the model that used *C*_{d} values calculated by the steady-flow CFD simulation (Fig. 6E), whereas the model with the *C*_{d} values from Hoerner (1965; Eqn 5 with *k*=0) decreased by 24% (to conform with the total drag force calculated by steady-flow CFD) and overestimated the resistive torque by 3.4% (Fig. 6F). In contrast, local torques along the length of the appendage yielded large differences between CFD for all three versions of the blade-element model (Fig. 6G–I). Nevertheless, a blade-element model of drag informed by a quantification of the total drag force in linear flow (Fig. 6F) offers an excellent approximation of the hydrodynamic resistance encountered by a rotating structure at the scale of a raptorial appendage. Such a model provided the basis of our calculations of the drag-torque index and the maximum torque and drag energy (see below).

### Drag measurements

Our drag measurements from dynamically scaled models allowed us to consider the hydrodynamics of different appendage shapes (Table S2). The two smasher taxa (*G. smithii* and *G. falcatus*) and the undifferentiated taxon (*H. californiensis*) exhibited higher drag when the dactyl was in the open position, while the drag on spearer appendages (*C. scolopendra* and *L. maculata*) was similar regardless of whether the dactyl was opened or closed against the propodus (Table S2, Figs 3, 5).

Our blade-element model succeeded in characterizing these drag measurements. The blade-element model calculated drag from measurements of the dimensions of the appendage (Eqn 5) and by fitting a value for the shape factor, *k*, included in the calculation of drag coefficient (Eqn 6). For each physical model, the blade-element model succeeded in characterizing the variation in drag with respect to flow velocity to generate a high coefficient of determination (*r*^{2}>0.97; Fig. 7). The smasher species, *G. smithii* and *G. falcatus*, with a closed dactyl, exhibited variation in drag with speed similar to that predicted for a uniform elliptical cylinder with the same mean dimensions. This result emerged in spite of the fact that these species exhibited shape factors (*k*≈0.10, Table S3) that were more than sevenfold greater than an elliptical cylinder (*k*=0.015). This discrepancy is explained by the substantial variation in cross-sectional shape along the length of a raptorial appendage. Conversely, the raptorial appendage of *C. scolopendra* was relatively uniform along its length (Fig. 3) and, consequently, was found to possess a shape factor (*k*=0.0153) that was similar to that of an elliptical cylinder.

### Effects of shape, size and kinematics on costs of movement

The drag-torque index provides a metric for the effects of the shape on the hydrodynamics of a raptorial appendage in rotation. Among the five species considered, we found that the two spearers had lower *T*_{d} than smashers and the undifferentiated *H. californiensis*. The spearers exhibited *T*_{d} values that were similar to that of a uniform elliptical cylinder (Table S3). This was because the spearers possessed a relatively uniform cross-sectional shape along the length of the appendage (Fig. 3). In contrast, the enlarged dactyl of a smasher at the distal end of the appendage yielded a relatively large value for the drag-torque index (Eqn 8). Opening the dactyl serves to move a portion of the appendage toward the distal end of the appendage, where drag generates a greater torque than when the dactyl is closed. As a consequence, the values of *T*_{d} in all species were greater with the dactyl opened than when it was closed.

We tested the individual effects of strike duration, dactyl length and drag-torque index over a series of 100 simulated strikes for three taxa (*G. smithii*, *L. maculata* and *C. scolopendra*) (Fig. 8). These simulations revealed that the maximum torque and drag energy were most sensitive to the size of the appendage and far less sensitive to strike duration and drag-torque index (Fig. 8).

## DISCUSSION

Drag generated by the strike of a mantis shrimp is dominated by size and velocity and less so by shape. Larger mantis shrimp, regardless of appendage type, move more slowly, which indicates the presence of trade-offs between angular kinematics and size during the diversification of mantis shrimp raptorial appendages. In spite of their distinct morphologies, differences in the shapes of the smasher and spearer appendages did not greatly alter drag. Nonetheless, we found that smashers reduced drag by closing the dactyl, and that spearers, even with very substantial dactyls, did not experience a large change in drag with the dactyl open.

A blade-element approach adequately approximated drag when compared with the more computationally intensive CFD approach. However, there were anomalies between the two approaches that warrant further investigation when considering more detailed, future analyses of energetic costs of raptorial strikes. In particular, the interaction between shape and span-wise flow suggests that shape may be important for local flow specific to rotational movement. Nonetheless, the reasonable correspondence between blade-element and CFD models is good news for coarse-grained future comparative analyses that are, from this point forward, more feasible than if CFD were required.

### Strike kinematics

The kinematic data from a medium-sized spearer, *C. scolopendra*, filled a key size range for the comparison of kinematics across stomatopod appendage types (Fig. 5, Table 1). Following the broad pattern across mantis shrimp, *C. scolopendra* struck with slightly slower kinematics than the smaller spearer *A. vicina*, and more slowly than the comparably sized smashers *G. smithii* and *N. bredini*. *Coronis scolopendra* used their spring and latch system, as indicated by the characteristic sliding and sudden outward rotation of the propodus and dactyl (Cox et al., 2014; deVries et al., 2012; Kagaya and Patek, 2016; Patek et al., 2004, 2007). Across species, angular kinematics decreased with increasing appendage size. Smashers exhibited faster kinematics for a given appendage size than non-smashers.

### Modeling approaches to the hydrodynamics of fast rotations

The simplified approach of a blade-element analysis adequately approximated hydrodynamic drag when compared with a computational fluid dynamic model. Our CFD model indicated regional variation in flow, pressure and hydrodynamic torque that was not predicted by the blade-element model (Fig. 6), which treats regions along the length of the structure as independently operating sources of drag. Despite these discrepancies, as long as the shape-dependent coefficients of drag used as input in the blade-element models do not overestimate drag in the case of steady linear flow, the blade-element model will also perform well to predict drag-torque for extremely fast rotations (Fig. 6E,F).

A previous study on the rotating head of a pipefish came to the same conclusion (Van Wassenbergh and Aerts, 2008), and our study now demonstrates that this approach is also valid for the more than 15-times higher angular velocities and more than 20-times higher angular accelerations of the fastest appendage strikes of mantis shrimp. A noteworthy discrepancy in the mantis shrimp modeling was that the pattern of local variation in the hydrodynamic torque along the length of the appendage was poorly predicted by the current blade-element models compared with the pipefish models (Fig. 6G–I versus fig. 9 in Van Wassenbergh and Aerts, 2008). Even the blade-element model of the mantis shrimp appendage that accounted for span-wise flow near the proximal and distal ends (under linear flow conditions) using locally reduced *C*_{d} values during the torque calculations still showed considerable local torque differences when compared with CFD (Fig. 6H). We assume this is caused by a different, yet related, simplification of blade-element models, namely the simplified treatment of local angles of attack.

Given that the current blade-element models treat the elements as elliptical cylinders, all surfaces in contact with the water are assumed to be parallel with the long axis. This is not the case in the more realistic CFD model that used frustum segments instead of elliptical cylinders. Therefore, the multi-frustum model appropriately accounted for variation in the angle of attack of the leading edge of the mantis shrimp appendage during rotation (Fig. 6C): the leading-edge surfaces at the distal end (frustum segments 16–18; Fig. 6A) move nearly parallel with respect to the surrounding water during rotation, whereas the more proximal leading-edge surfaces (frustum segments 10–15; Fig. 6A) move nearly perpendicularly to the water. This is reflected in the pressure patterns on the surfaces: no positive pressures are exerted by the water on the distal part of the leading edge, whereas a large zone of positive pressure is located on the more proximal part of the leading edge.

Consequently, the larger discrepancy in local torques in the mantis shrimp model (Fig. 6A) compared with the pipefish is probably due to the larger variation in the angles of attack at the distal part of the rotating structure. This effect will also be less prominent in the spearer mantis shrimp, because of their lower variation in chord lengths along their appendages (Fig. 3). Despite these local torque discrepancies, the overall congruence between the total torque output from CFD and blade-element modeling is significant for future comparative studies because it demonstrates that hydrodynamic effects that are inherent to rotational flow do not substantially impact the overall drag-torque over a wide range of angular velocities and accelerations.

Both CFD and blade-element approaches present limitations on how we have modeled drag. These techniques do not account for the effects of a transition to turbulent flow in the boundary layer of the structure. This transition, which for smooth surfaces in steady flow occurs at a *Re* of approximately 200,000, would cause drag coefficients to decrease significantly (Hoerner, 1965). However, this critical *Re* was exceeded by one species (*G. smithii*) and only for a brief period upon attaining maximum speed at the distal tip of the appendage (Table 1). A turbulent transition even in this case appears unlikely. The flow over the appendage is predominantly chord-wise, which suggests that the chord is the appropriate characteristic length for the *Re* that is predictive of the turbulent transition, which would likely yield a *Re* value below 200,000.

Our models additionally neglect surface roughness, which may induce a turbulent transition in steady flow at *Re*>20,000 (Hoerner, 1965). However, flows are far from steady during the strike of mantis shrimp. Studies on flow accelerated from rest in pipes show that the transition to a turbulent boundary layer is delayed and thereby causes the transition to occur at significantly higher *Re* (Lefebvre and White, 1989). Comparing the wake of our accelerated simulations with steady flow simulations from CFD support the notion that the flow corresponds to the patterns of a lower *Re* (Fig. S2). Unfortunately, the combined effect of surface roughness (decreasing critical *Re*) and acceleration from rest (increasing critical *Re*) on the transition to a turbulent boundary layer is unknown. Because accelerations are extremely high, and the appendage surfaces relatively smooth, we have assumed that the boundary layer remains laminar for mantis shrimp strikes.

The suite of approaches taken in this study, including our drag measurements, blade-element analyses of shape-variant mantis shrimp appendages, and CFD analyses of rotating multi-frustum models, lends strong support to a simplified approach going forward. Furthermore, the similar excursions of strikes across taxa (Fig. S1) allow an additional level of comparative analysis that incorporates a similar kinematic profile and thus reduces the need for time-intensive kinematic analyses. Future comparative studies could therefore be performed based primarily on morphology and thereby build on and integrate the substantial understanding of morphological and mechanical evolution of this system (Anderson et al., 2014; Anderson and Patek, 2015; Blanco and Patek, 2014; Claverie et al., 2011; Claverie and Patek, 2013; Patek et al., 2007, 2013; Rosario and Patek, 2015).

### The evolution of the raptorial appendage

Shape variation did not play a major role in drag, which suggests that the notable differences in smashing and spearing appendages can be largely attributed to their divergence in predatory function and morphological robustness for impaling and crushing prey (Anderson et al., 2014, 2016b; Claverie et al., 2011; Claverie and Patek, 2013; deVries et al., 2012). In only one dimension did shape play a notable role: regardless of whether the dactyl was open or closed, the spearers had lower torque resulting from drag (*T*_{d}) than smashers. In contrast, the smashers experienced greater drag energy and drag torque with an open dactyl than with the dactyl closed in a hammering position.

The combined effects of size and speed on drag generation were substantial across mantis shrimp. A spearing mantis shrimp's strike encounters drag energy and torque comparable to that of a smashing mantis shrimp that is an order of magnitude smaller (Fig. 8B,C). If one considers drag as a proxy for energetic costs, then an upper limit on energy costs may explain the apparent evolutionary limits on maximal body size in smashers (Blanco and Patek, 2014; Patek, 2015). In other words, to move extremely fast, smasher appendages must be small. The energetic costs may be even more stringent than evidenced in the present paper, which allowed mathematical models of larger appendages to rotate more quickly; a previous hydrodynamic analysis demonstrated that higher displacement strikes intuitively associated with greater speeds (given the longer out-lever of the appendage) actually move more slowly than expected, because of the energetic costs of moving through water (McHenry et al., 2012). Future comparative analyses that more fully account for energetic costs may pinpoint more subtle variation within and across species that balances the trade-offs among these key variables.

### Conclusions

Unlike the hydrodynamics of locomotor systems that must balance propulsion with drag reduction, the diversification of feeding structures operate under a different constellation of factors, such as fracture resistance, effective prey capture morphology and the reduction of feeding costs (Anderson et al., 2016a; Anker et al., 2006; Full et al., 1989; Vermeij, 1987; Weaver et al., 2012). This study demonstrates that the shapes of fast, centimeter-scale predatory structures can diversify with relatively minimal costs within the steep constraints of size and kinematics in the generation of drag.

## Acknowledgements

We thank E. Murphy and P. Green for their help collecting specimens, and two anonymous reviewers for their constructive feedback. We are grateful for the hospitality and advice offered by R. Heard during our fieldwork with *C. scolopendra*. Micro-CT scans were performed at the Center for Nanoscale Systems (CNS) at Harvard University, a member of the National Nanotechnology Infrastructure Network (NNIN) (NSF ECS-0335765).

## Footnotes

**Author contributions**

M.J.M. contributed to the design of the project, performed mathematical modeling analyses, made figures, and contributed to text generation and editing. P.S.L.A. collected drag flume data, performed drag and PGLS analyses, made figures and contributed to writing sections of the paper. S.V.W. performed the computational fluid dynamics analyses and contributed to figures and text. D.G.M. collected and analyzed *C. scolopendra* kinematics and contributed to figures and text. A.S. helped collect drag flume data and print models. S.N.P. conceived of the project, performed the cross-species kinematic analyses, made figures and wrote the paper.

**Funding**

Funding was provided by National Science Foundation grants to S.N.P. (IOS-14391050), M.J.M. (IOS-13541042 and IOS-0952344), A.S. (IOS-1256602) and the John Simon Guggenheim Memorial Foundation to S.N.P. S.V.W. was supported by the Agence Nationale de la Recherche (ANR-16-ACHN-0006-01).

**Data availability**

All datasets and computer coding are available from the Dryad Digital Repository http://dx.doi.org/10.5061/dryad.578j2

## References

**Competing interests**

The authors declare no competing or financial interests.