Significance Plunge-diving is a very unique foraging method in the animal kingdom. A limited set of water birds exhibit this behavior, and only one family of seabirds (Sulidae) exhibit this behavior at high speeds. We studied the stability of the bird’s slender and seemingly fragile neck during a plunge-dive by conducting simple experiments that mimic this behavior. An elegant analysis of the interaction among hydrodynamic forces, neck elasticity, and muscle contraction reveals that seabirds dive at appropriate speeds to avoid injury. Considering the popular recreational sport of diving, we also find a diving speed limit for humans to avoid injury.

Abstract In nature, several seabirds (e.g., gannets and boobies) dive into water at up to 24 m/s as a hunting mechanism; furthermore, gannets and boobies have a slender neck, which is potentially the weakest part of the body under compression during high-speed impact. In this work, we investigate the stability of the bird’s neck during plunge-diving by understanding the interaction between the fluid forces acting on the head and the flexibility of the neck. First, we use a salvaged bird to identify plunge-diving phases. Anatomical features of the skull and neck were acquired to quantify the effect of beak geometry and neck musculature on the stability during a plunge-dive. Second, physical experiments using an elastic beam as a model for the neck attached to a skull-like cone revealed the limits for the stability of the neck during the bird’s dive as a function of impact velocity and geometric factors. We find that the neck length, neck muscles, and diving speed of the bird predominantly reduce the likelihood of injury during the plunge-dive. Finally, we use our results to discuss maximum diving speeds for humans to avoid injury.

Nature contains several species of creatures that interact with the air–water interface (1). A number of bird species are able to dive into water from the air as a hunting mechanism (e.g., kingfishers, terns, and gannets), a behavior known as plunge-diving (2, 3). Some seabirds, like the northern gannet, are highly specialized plunge-divers, making 20–100 dives per foraging trip, diving from heights of 5–45 m, and attaining speeds of more than 20 m/s (4⇓⇓–7). Thus, the bird’s structure and behavior have presumably evolved to withstand a variety of high dynamic stresses, because no injuries have been reported in plunge-diving seabirds. Biologists have previously focused on the diving behavior in terms of ecological factors, such as diving depths, prey species, and hunting success rate (8⇓–10), and physiological features, such as the role of vision while crossing the air–water interface (11, 12). Unique kinematic and morphological features during the dive have also been observed, such as having a sharp, arrow-like body posture and a straight, long, and slender neck (13, 14). However, a mechanical understanding of plunge-diving birds is not well-established.

To study such a phenomenon, Morus bassanus (hereafter gannets) and Sula leucogaster (hereafter boobies), from the Sulidae family, are used as a model species due to their highly specialized diving characteristics (5, 13). First, they plunge-dive at very high speeds, using that momentum to carry them to some depth. Then, they use their webbed feet and/or wings to propel themselves further underwater, like penguins and cormorants (15, 16). Although plunge-diving at high speeds allows the bird to dive deeper, it induces much larger stresses on the seabird’s body than pursuit diving alone (13). The two main forms of plunge-diving observed are known as the V-shaped dive and the U-shaped dive (5). During V-shaped dives, the seabird impacts the surface at an angle, whereas during U-shaped dives the impact trajectory is more perpendicular to the surface (14, 17). Although the mechanical forces may differ between the two dives, both U-shaped and V-shaped dives experience an axial force significantly larger than a transverse force. Therefore, this present study focuses on the U-shaped dive, which is a model for understanding the effect of an axial force on the risk of a buckling neck.

From a mechanics standpoint, an axial force acting on a slender body may lead to mechanical failure on the body, otherwise known as buckling. Therefore, under compressive loads, the neck is potentially the weakest part of the northern gannet due to its long and slender geometry. Still, northern gannets impact the water at up to 24 m/s without injuries (18) (see SI Appendix, Table S1 for estimated speeds). The only reported injuries from plunge-diving occur from bird-on-bird collisions (19). However, for humans, diving into water at speeds greater than 26 m/s risks severe fractures in the cervical or thoracic vertebrae and speeds greater than 30 m/s risk death, regardless of impact orientation (20⇓⇓⇓⇓⇓–26). Understanding the bird plunge-dive may further explain methods of injury prevention in human diving.

In this present study, we investigate how birds are able to dive at high speeds and sustain no injury, given the morphology of the head and neck. Due to its long, slender geometry, the seabird’s neck is the region with the greatest potential for mechanical failure or instability under high dynamic loading. In reduced-order experiments, we simplify the seabird system as a long, thin, elastic beam attached to a rigid cone, which represent the bird’s neck and head, respectively. By modeling the bird’s neck as an elastic beam, we can use the buckling and nonbuckling behaviors of the elastic beam to represent the stability of the seabird’s neck. A linear stability analysis is used to obtain a theoretical prediction of the buckling transition. The effect of neck muscles is discussed in terms of modifying the buckling criterion. We then show that plunge-diving seabirds have a unique morphology, appropriate diving speeds, and strong neck muscles that will allow them to dive safely at high speeds.

Results Plunge-Diving Seabirds. To characterize the plunge-diving mechanism of seabirds, a salvaged northern gannet is prepared in the diving posture and is released into a water tank as shown in Fig. 1A (Materials and Methods). Upon water entry, in which momentum carries the bird through the water (11), three different phases become apparent: (i) the impact phase, (ii) the air cavity phase, and (iii) the submerged phase, which is characteristic of a classic water entry problem (27). The impact phase occurs when the tip of the beak first makes contact with the water surface until the head becomes submerged ( t < H head / V ; H head is the head length and V is the diving speed at impact). The air cavity phase occurs between the head’s being submerged until the seabird’s chest makes contact with the water surface [ t < ( H head + L neck ) / V ; L neck is the neck length, which is also close H head ]. The submerged phase occurs after the bird’s chest impacts the water, closing the air cavity. Fig. 1. (A) A deceased, frozen northern gannet impacts the water surface vertically at V ≃ 5.5 m/s and develops an air cavity around its neck as it descends (Movie S1). The impact phase occurs in the range when the tip of the beak first makes contact with the water surface until the head becomes submerged. The air cavity phase occurs between the head’s being submerged up until the seabird’s chest makes contact with the water surface. The submerged phase occurs when the bird’s chest impacts the water, closing the air cavity. Note that the head length ( H head ) and the neck length ( L neck ) are approximately the same. (B) Top and side views of a northern gannet (M. bassanus) skull and a brown booby (S. leucogaster) skull. Arrows indicate the location of the naso-frontal hinge. The air cavity phase [ H head / V < t < ( H head + L neck ) / V ] is the most interesting because it provides the greatest potential for neck injuries during the plunge-dive. During the air cavity phase, the head is subjected to hydrodynamic forces causing the head to decelerate while the rest of the body continues to descend downward. This results in an axial compressive load on the neck. Once the chest makes contact with the water in the submerged phase, the compressive load on the neck will be alleviated due to the hydrodynamic force on the chest. A nondimensional number representing a ratio of hydrodynamic drag to the neck’s elasticity [ η = ρ f V 2 R head 2 / ( E I L neck − 2 ) ; ρ f is the fluid density, R head is the head radius, E is the elastic modulus of the overall neck, and I is the area moment of inertia of the neck; Materials and Methods and SI Appendix, Fig. S1] is about 4.7 for birds diving at a speed of 24 m/s. This simple scaling analysis shows that the drag force may exceed the compliance of the neck, potentially leading to injuries. However, for a better assessment of the neck stability, we need to examine the effect of head shape and neck muscles. Physical Experiment to Mimic Plunge-Diving. To further explore this fluid–neck interaction, we design a reduced-order experiment by approximating the neck (a composition of bone, muscle, and skin) as an elastic beam and the head as a rigid cone. A cone–beam system was fabricated to effectively model the head–neck interaction during impact (Materials and Methods). Various geometric parameters (i.e., cone angle, cone radius, and beam length) and impact velocities were tested, producing a range of drag to elasticity ratio to be η = O ( 10 − 2 − 10 2 ) , which encompasses the drag to elasticity ratio value for plunge-diving birds. Fig. 2A shows the dynamics of a cone–beam specimen remaining stable through the air cavity phase, H cone / V < t < 2 H cone / V . Fig. 2B shows a specimen with a longer beam length becoming unstable through the air cavity phase. The effect of other parameters (cone angle, velocity, and length) on the stability of the beam can be seen in SI Appendix, Fig. S2. In general, an increase in velocity and beam length will increase the likelihood of buckling. Fig. 2. Time sequence images of cone–beam specimens impacting the water surface as inspired by a diving seabird. Here, β = 30°, V = 0.65 m/s, and the moment of impact is set at T ∼ = 0 . (A) A specimen with beam length L = 5 cm exhibits a stable, nonbuckling behavior (Movie S2). (B) A specimen with beam length L = 8 cm exhibits unstable, buckling behavior. (C) Nondimensional change in amplitude ( Δ Y / h ) vs. nondimensional time ( T ∼ ) for cases A and B. Case A does not exceed Δ Y / h = 1 (stable), whereas case B exceeds Δ Y / h = 1 (unstable). (Inset) Increasing velocity will increase the buckling amplitude. Forces on Cone or Head. Hydrodynamic drag is the main force that acts on the cone during impact. When the cone first enters the water (during the impact phase), the drag force is primarily induced by a change in added mass. Hence, the drag force increases in time during the impact phase; F Drag ( t < H cone / V ) = 2 ρ f V 4 tan 3 ( β ) t 2 . During the air cavity phase, the drag force and hydrostatic pressure force can be expressed as F Drag ( t ≥ H cone / V ) = π 2 ρ f C d R cone 2 V 2 tanh ( β ) , and F Hydr ( t ) = π ρ g R cone 2 ( z ( t ) − 2 3 H c o n e ) , where β is the cone half-angle, R cone is the cone radius, and z ( t ) is the distance between the cone tip and the free surface (Materials and Methods). Here, the drag coefficient, C d , is chosen to be 2/3 as a fitting parameter. Smaller cone angles will reduce the hydrodynamic drag forces during the impact phase, but during the air cavity phase the hydrostatic pressure will increase more rapidly due to the larger cone height. Force data on a cone with β = 12.5 ° are presented in Fig. 3A. Measured force and time are normalized as F ∼ = F / ( π 2 ρ f R cone 2 V 2 tan β ) and T ∼ = t / ( H cone / V ) , respectively. The measured force initially rises parabolically during the impact phase due to the strong time dependence in F Drag ( t < H cone / V ) . At later times ( t ≥ H cone / V ), the drag force becomes constant. However, the hydrostatic pressure force linearly increases during the air cavity phase, which is predicated by the analytical models discussed above. So, when the cone has a larger cone height and impacts at a low speed, the hydrostatic pressure force plays a larger role because more time is needed for the cone to reach T ∼ = 2 . Shorter cone heights with higher speeds are less affected by F Hydr at T ∼ = 2 . Fig. 3. Force on the cone and 3D-printed northern gannet skull during impact. Force data are collected at four different impact velocities ranging from 2.1 to 3.2 m/s. The forces are normalized by ∼ ρ f R 2 V 2 . (A) (Upper) Time sequence of a cone of R = 3 cm and β = 12.5 ° entering water at 2.4 m/s. (Lower) Nondimensionalized force and time of experimental data for the cone cases. (B) (Upper) Time sequence of a 3D-printed northern gannet skull entering the water at 3.1 m/s. (Lower) Nondimensionalized force and time of experimental data for the skull cases (Movie S1). Next, we consider the force on a 3D-printed skull of a northern gannet. Based on the geometry of the skull, three distinct sections are identified (Fig. 1B and SI Appendix, Fig. S3). The first section is between the tip of the beak to its base, where a hinge [naso-frontal hinge (28)] runs along the dorso between the beak and the forehead (Fig. 1); the second section is between the naso-frontal hinge and small protrusions near the back of the skull (zygomatic process of the Os squamosum); the third section is between the protrusions and the end of the skull (Prominentia cerebellaris) (29). Assuming that the skull is two cones of different angles in tandem, the force measurement during the impact phase shows two distinct time-dependent curves as predicted by our analytical expression described above (Fig. 3B). This result indicates that the axial force acting on the neck of the plunge-diving bird increases with the skull radius, the impact velocity, and, most importantly, the beak angle. Transition to Buckling. The transition from stable to unstable beams depends on the impact velocity, geometric factors, and material properties of the beam and the cone. The critical compressive force to buckle the beam is calculated from a linear stability analysis resulting in the dispersion relation. In order for the beam to buckle, the highest growth rate at some given time must lie in the unstable region (Fig. 4A), in our case nondimensional wavenumber greater than π ( k L ≥ π ) (30, 31). In other words, buckling only occurs when the most unstable wavelength is shorter than the beam length. At the moment of the fully submerged neck ( t = 2 H head / V , or T ∼ = 2 ), we obtain a buckling criterion based on the spatial condition for the beam as 2 A c E ρ b ρ f [ F Bend − F Hydr + F W ] < V c C d tanh ( β ) , [1]where V is impact velocity, β is cone half-angle, c = E / ρ b is the speed of sound in the material, E is beam elastic modulus, ρ b is beam density, ρ f is fluid density, A c is projected cone area, C d is drag coefficient, F Bend is bending force, F Hydr is hydrostatic pressure force, and F W is cone weight (SI Appendix, Fig. S4). From a biological perspective, choosing to analyze the beam behavior at T ∼ = 2 is analogous to the time when the bird’s chest would impact the water, or when the compressive load on the neck begins to be alleviated. Fig. 4B shows the region of stability, which predicts the stability of the beam under various conditions, and is in good agreement with experimental data. Fig. 4. (A) Growth rate vs. nondimensionalized wavelength. Each curve represents a different time. The black curve is the moment when t = H cone / V . The red curve is the moment when t = 2 H cone / V . (B) Phase diagram of the cone–beam system. Various shapes represent different cone half-angles. Blue markers indicate stable, nonbuckling cases; red markers are unstable, buckling cases; light-blue markers represent the specimens in the transition regime, exhibiting buckling and nonbuckling behaviors under the same test conditions. Both the brown booby and northern gannet plunge-dive in the stable region of the phase diagram. (C) CT scan of the northern gannet. A bundle of muscles is circled in white behind the skull. (D) Contracting muscles help to keep the neck straight and therefore act as a stabilizing mechanism during the plunge-dive. Using morphological and material properties obtained from the salvaged bird, we find that the plunge-diving birds dive in the stable region of the transition diagram. However, this analysis neglects the effects of the neck muscles, which leads to another question addressed in the next section. What role does the neck musculature play in preventing neck injuries during the plunge-diving behavior? Effect of Muscles. The motion and strength of an animal’s neck result from the coupling between bone and muscles (32, 33). The total force generated by the muscle bundle will depend on the length and cross-sectional area of muscle fibers. Neck muscles in plunge-diving birds are mostly concentrated near the head and the thorax of the bird, as shown in Fig. 4C. The muscles connect the body, the vertebrae, and the skull by a series of thin muscle sheets and tendons (34). Additionally, the necks of gannets and boobies, similar to those of other birds, have an S shape, due to their vertebrae morphology and connecting design (34, 35), increasing the complexity of biomechanical analysis. One may note that the S shape would serve as a primary mode of buckling in the neck. However, the fact is that the musculature plays an important role in stabilizing a straight neck as a whole, and also maintaining the S shape of the spine. Therefore, we simplified the complex network of muscles using segmentation and reconstruction of computed tomography (CT) images of the lateral, dorsal, and ventral musculature (Fig. 4C). By muscle contraction, the tendons put some stabilizing tension on the bones, straightening the neck out and fixing the bones into place before the impact. We approximate the effect of the muscles as a continuously distributed load acting tangentially along the neck (Fig. 4D) (36). The muscle force f muscle per cross-sectional area is estimated by measuring the cross-sectional area of the neck muscles in CT-scanned images (37). Then, we find that a resisting bending force is modified as F Bend = ( 2 π 2 E I ) 1 / 3 f muscle 2 / 3 due to the contracting muscle force. In our stability analysis, the inclusion of the muscle force will place the plunge-diving birds higher in the transition diagram, thus further away from the buckling transition line. The critical impact force to overcome this modified muscle/bending force was estimated to be 3.4 kN, which is two orders of magnitude higher than the 30 N produced by the combined hydrostatic and drag force, thus allowing the bird to dive safely at high speeds.

Discussion The results help to reveal the mechanisms (in addition to visual accommodations) by which plunge-diving birds are able to dive at incredibly high speeds with no injuries (19). This is primarily attributed to the neck length and chosen diving speeds, which stay in parameter regimes that prevent the neck from bending under compressive loads. The neck muscles move plunge-divers further away from the buckling transition. In fact, it would take about 80 m/s for the plunge-diving seabird to sustain a neck injury based on our analysis. Furthermore, this study may elucidate safe diving speeds for humans. We consider feet-first dives, which gives a higher survival rate (20). Human feet are flat with large surface areas; average foot areas for males and females are 0.06 and 0.05 m2 (38), respectively. At a diving speed of 24 m/s, the compressive force that a human would experience is about 14 kN, which well exceeds a range of maximum compressive forces (0.3–17 kN) to cause neck injury (39). The impact force exceeds the critical maximum compressive force (17 kN) at a diving speed of about 24 m/s (for trained individuals, i.e., stunt divers). This critical diving speed is consistent with the maximum speed (≈26 m/s) for spinal fractures reported in case studies (SI Appendix, Table S2).

Materials and Methods Salvaged Bird. A salvaged northern gannet (M. bassanus) was obtained for analysis. The elastic modulus of the neck ( E ≃ 8.6 MPa) was determined based on the neck’s curvature from an applied load (SI Appendix, Fig. S1B). The bird was then frozen in a position so that the neck was extended straight. To understand morphological properties, the frozen bird was CT-scanned (Toshiba Aquilion 16; 100 kVp, 125 mA, 0.5-mm slice thickness, 512 × 512 reconstruction matrix) and the resulting images provided us with the neck length ( L neck ≃ 21 cm), neck and head radius ( R neck ≃ 2 cm and R h e a d ≃ 2.5 cm), and skull angle (β = 11°) (reconstruction and visualization of images were done using Horos open-source software, v. 1.0.7, https://www.horosproject.org/). Therefore, the bending rigidity becomes E I ≃ 1.1 N ⋅ m2. Segmentation and reconstruction of the skeleton and musculature of the neck were done using Mimics software (Materialise NV), providing the cross-section areas for the muscle force calculations. The frozen bird was then subjected to a series of drop tests into a tank of water. High-speed footages of the impact, air cavity, and submerged phases were acquired. Skull Specimens. Several skull specimens for different species of gannets and boobies were acquired from the Smithsonian Museum of Natural History (M. bassanus, n = 14; Morus capensis, n = 5; Morus serrator, n = 2; Sula dactylatra, n = 2; Sula sula, n = 3; and S. leucogaster, n = 3). Two distinct regions are noted: (i) between the tip of the beak to the naso-frontal hinge having a half-angle of β 1 = 7.9° ± 0.6° and (ii) between the naso-frontal hinge to the zygomatic process (from the squamous bone) with β 2 = 12.3° ± 0.8°. The average skull radius was measured to be 2.4 ± 0.3 cm. Details of measurements can be seen in SI Appendix, Fig. S3. Muscle Cross-Sectional Area. The northern gannet’s neck musculature was divided in dorsal, ventral, and lateral (Fig. 4C; red for dorsal muscle and yellow for ventral muscle). Mesh masks were created using threshold selection (Hounsfield unit: −140, 299) and cleaned for segmentation and reconstruction in Mimics program. To measure the cross-sectional area of the dorsal and ventral musculature, the neck was divided in anterior (five vertebrae close to the skull) and posterior (vertebrae 9–14, closer to the thorax) portions; vertebrae 6–8 comprise the midsection of the S-shaped neck. After segmentation, seven points along the anterior and posterior portions of the neck were selected to extract the cross-sectional area measurements. Mask area measurements (square millimeters) of 10 sequential slice images were averaged for each of the seven points. Values from the seven points were averaged (anterior musculature: ventral, 195.87 ± 41.22 mm2 and dorsal, 319.21 ± 186.78 mm2; posterior musculature: ventral, 88.79 ± 23.76 mm2 and dorsal, 181.40 ± 68.44 mm2) to determine musculature forces to avoid neck buckling. Physical Experiment. To simulate the plunge-diving seabird’s head–neck coupling, a cone–beam system was developed (SI Appendix, Fig. S5). Cones with radii of 1.27 cm and 2.0 cm were either 3D-printed (Makerbot Replicator 2X, ABS plastic) or manufactured (acrylic). The cone half-angles, β, were 12.5°, 30°, 35°, 40°, 45°, 50°, and 58°. Rectangular elastic beams were created using vinylpolysiloxane (Elite Double 22; Zhermack Co.) (E = 0.95 MPa and ρ b = 1,160 kg/m3). Whereas a bird’s neck has a circular cross-section, the elastic beams consist of a rectangular cross-section to control the beam’s bending plane for image analysis. The elastic beam is attached to the cone on one end and clamped at the other end at some distance L (2–10 cm) from the cone base. Using the MATLAB image processing toolbox, the amplitude of the beam ( Δ Y) was correlated with the change in distance between the cone and the clamp ( Δ Z). The cone–beam system is dropped from various heights, resulting in impact velocities ranging from 0.5 to 2.5 m/s, and recorded using a high-speed camera (IDT-N3, 1,000 frames per s). At least five trials are conducted for each set of the experimental parameters. The changing vertical distance ( Δ Z) is calculated during a time frame ranging from slightly before impact to a time when the specimen is submerged a distance equivalent to two cone heights below the water surface. By measuring Δ Z, we can determine the amplitude Δ Y while avoiding interference of the water splash. The amplitude is nondimensionalized using the beam’s thickness, Δ Y / h , from which we can sort the data into three separate categories of stability: buckling, transition, and nonbuckling. After processing all high-speed videos for both amplitude and velocity data, our experiments exhibit three states: stable, unstable, and transitionally unstable. Quantitatively, these states can be characterized by distinct ranges of the nondimensional amplitude. The stable state is characterized by a nondimensional amplitude range less than one, which corresponds to the nonbuckling behavior of the beam; conversely, the unstable state has a nondimensional amplitude greater than one, which corresponds to the unstable buckling behavior. After repeated trials of a single case, the case is considered stable if fewer than 20% of the trials buckle. If more than 80% of the trials buckle, then the case is unstable. If 20–80% of the trials buckle, then the case is characterized as transitionally unstable. Derivation and Measurement of Forces. The drag force during the impact phase is derived from the Euler–Lagrange equation. The Lagrangian is defined as ℒ = K . E . − P . E . The kinetic energy term is described as K . E . = 1 / 2 ( m cone + m add ) V 2 , and the potential energy term, P . E . , is neglected because it is small compared with the kinetic energy term. The added mass term is dependent on the instantaneous radius of wetted area, m add = 4 / 3 ρ f r ( t ) 3 . Further relationships to consider are r ( t ) = z ( t ) tan ( β ) and z ( t ) = V t . Now, the Euler–Lagrange equation, ∂ ℒ / ∂ z − d / d t ( ∂ ℒ / ∂ z ˙ ) = 0 , can be reduced to F Drag ( t < H cone / V ) = 2 ρ f V 4 tan 3 ( β ) t 2 . After the impact phase, the drag force is no longer time-dependent and simply becomes F Drag ( t ≥ H cone / V ) = π / 2 ρ f C d R 2 V 2 tanh ( β ) . The hydrostatic pressure force was determined by integrating the pressure along the surface of the cone, resulting in F Hydr ( t ) = π ρ g R cone 2 ( z ( t ) − 2 / 3 H c o n e ) . A rigid steel rod connects the cone/bird skull to the force transducer (LCM-105-10; Omegadyne, Inc.). The force transducer is connected to a signal conditioner (2310; Vishay), which collects data at a sampling rate of 1 kHz. The high-speed camera was used again to determine the impact velocity at 1,000 frames per s. At least five trials were taken for the cone with β =12.5° and the 3D-printed northern gannet skull impacting the water from 2.0 to 3.2 m/s (SI Appendix, Fig. S6). Derivation of Dispersion Relation. Under an axial force, the lateral displacement, Y ( Z , t ) , of a slender elastic beam can be described by the linearized Euler–Bernoulli beam equation, ρ b A b ∂ 2 Y ∂ t 2 + ∂ ∂ Z ( F ( t ) ∂ Y ∂ Z ) + E I ∂ 4 Y ∂ Z 4 = 0 , where ρ b is the beam density, A b is the cross-sectional area of the beam, E is the elastic modulus, I is the area moment of inertia of the beam, and F ( t ) is the axial load. By assuming the normal mode [ Y ( Z , t ) = Y 0 e ω t + i k Z ] of the beam deflection, the dispersion relation was determined to be ω 2 = E I ρ b A b k 2 ( F ( t ) E I − k 2 ) . It is noteworthy that the axial force [ F ( t ) ] depends on time and therefore the growth rate (ω) of perturbations changes over time. Assuming that different wavemodes are independent, the growth rate, ω, at any given time is dependent only on the history of the growth rate. Therefore, the growth rate with a time-dependent force can be described by integrating ω over time at a given wavelength (k). A similar approach has been used previously in the Rayleigh–Plateau instability of a crown splash (30, 31). We use the time-dependent force, which has a discontinuity at T ∼ = 1. Fig. 4A shows the dispersion relation between the integrated growth rate ( ∫ ω d t ) and nondimensional wavelength ( k L ). At any instance, the most unstable mode ( k L ) is determined by finding the maximum of the integrated growth rate [ ∂ ( ∫ ω 2 d t ) / ∂ k = 0 ], which is marked in solid circles. We find that the most unstable mode increases in the beginning until T ∼ = 1. This trend of the increasing unstable mode is anticipated primarily due to the increasing compressive force due to impact. Beyond T ∼ = 1, the most unstable mode either shifts to a lower k L because steady drag is greater than the impact force (SI Appendix, Fig. S7A; see β = 50°) or continues to increase because the hydrostatic pressure force is greater than the steady drag (SI Appendix, Fig. S7B). Spatial Stability. If we assume force is time-independent (SI Appendix, Fig. S8), the most unstable mode is given by ∂ ω 2 / ∂ k = 0 . Here, we find the critical wavemode to be k crit = F / ( 2 E I ) . When the most unstable wavelength is less than the beam length ( 2 π / k crit < L ), the beam is unstable (buckled). As a result, the beam will be unstable during impact when 8 E I π 2 / F ( t > H / V ) < ( α L ) 2 , where α = 2 is the prescribed boundary condition in which one end of the beam is fixed and the other end is free to move laterally. After substituting the forces at the moment of T ∼ = 2, the resulting equations yield F Drag + F Hydr − F W > F Bend ( = 2 π 2 E I L 2 ) , [2] where F Bend , F Hydr , and F W are the resistive bending force, hydrostatic pressure force, and weight of the cone, respectively. Rearranging terms in the above equation yields Eq. 1. In experiments, we chose a reference moment for the instability of the beam as when the cone reaches two cone heights (or t = 2 H / V ) below the free surface. Therefore, we incorporate the impact force into our analysis, which produces an instability criterion in terms of geometric factors, material properties, and impact velocity. Neck Muscle Resistance. The results from the previous sections indicate that the plunge-diving seabirds are able to dive safely at high impact speeds. The analysis, however, neglects the role of neck muscles. We consider the neck muscle force as a distributed follower load acting tangentially along the beam: F = F Drag + F Hydr − F W − ∫ f muscle d Z . [3] To simplify the analysis, the muscle force can be approximated as a constant force per unit length, f muscle . The muscle force per cross-sectional area is estimated based on the value in ref. 37: 37 N/mm2 for geospiza fortis and 17 N/mm2 for geospiza fuliginosa. In this study, we chose the lowest value (17 N/mm2) available in the literature to estimate the force generated by the anterior neck muscle as f muscle ≃ 17 N/mm2 × (196 + 319) mm2/ L neck ≃ 4.2 × 104 N/m. There is a critical length for which a beam retains neutral stability when an axial force competes with a distributed follower load (36). In our case, this length is defined by L N = ( F Drag + F Hydr − F W ) / f muscle . From our linear stability analysis, we have F Drag + F Hydr − F W = 2 π 2 E I / L N 2 , which will yield F Drag + F Hydr − F W = F Critical = ( 2 π 2 E I ) 1 / 3 f muscle 2 / 3 . If the hydrostatic and drag force exceeds the critical force, then the neck muscles will not be able to withstand the hydrodynamic forces, causing the bird neck injuries during the plunge-dive. By including the effects of the neck muscles, we speculate that the plunge-diving birds will move even further from the transition line in the stability diagram.

Acknowledgments We thank Alex Ochs, Grace Ma, Andrew Marino, Thomas Moore, and Yuan-nan Young for their initial contributions. This work was partially supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico Grant 246819/2013-8 (to L.S.), Virginia Tech Institute for Critical Technology and Applied Science, and National Science Foundation Grants CBET-1336038 (to B.C., S.G., and S.J.) and PHYS-1205642 (to S.G. and S.J.).