Penguins are adapted to underwater life and have excellent swimming abilities. Although previous motion analyses revealed their basic swimming characteristics, the details of the 3D wing kinematics, wing deformation and thrust generation mechanism of penguins are still largely unknown. In this study, we recorded the forward and horizontal swimming of gentoo penguins (Pygoscelis papua) at an aquarium with multiple underwater action cameras and then performed a 3D motion analysis. We also conducted a series of water tunnel experiments with a 3D printed rigid wing to obtain lift and drag coefficients in the gliding configuration. Using these coefficients, the thrust force during flapping was calculated in a quasi-steady manner, where the following two wing models were considered: (1) an ‘original’ wing model reconstructed from 3D motion analysis including bending deformation and (2) a ‘flat’ wing model obtained by flattening the original wing model. The resultant body trajectory showed that the penguin accelerated forward during both upstroke and downstroke. The motion analysis of the two wing models revealed that considerable bending occurred in the original wing, which reduced its angle of attack during the upstroke in particular. Consequently, the calculated stroke-averaged thrust was larger for the original wing than for the flat wing during the upstroke. In addition, the propulsive efficiency for the original wing was estimated to be 1.8 times higher than that for the flat wing. Our results unveil a detailed mechanism of lift-based propulsion in penguins and underscore the importance of wing bending.
Among the wing-propelled diving birds, penguins (Spheniscidae) are considered to be the most specialized for underwater swimming because they are free from the constraints of aerial flight (Elliott et al., 2013; Storer, 1960). Their excellent swimming ability includes high-speed foraging, long migration and deep diving. For example, field studies have shown that the average foraging speed of gentoo penguins (Pygoscelis papua) is 2.3 m s−1 (Sato et al., 2007), that Fiordland penguins (Eudyptes pachyrhynchus) swim 80 km per day (Mattern et al., 2018) and that the diving depth of emperor penguins (Aptenodytes forsteri) reaches 564 m (Wienecke et al., 2007). Such excellent performance is realized by flapping-wing propulsion adapted to the physical properties of water. Since the density of water is more than 800 times greater than that of air, the buoyancy on the body and hydrodynamic loads on the wing in penguins are much greater than those in aerial birds. Penguins have to generate a vertical downward force to counteract the upward buoyancy (Sato et al., 2002) and their wings have thick aerofoil sections (maximum thickness: 17.5% for gentoo penguin) supported by arm and hand bones (Bannasch, 1995). In contrast, aerial birds generate a vertical upward force mainly by the downstroke to support their body weight in the air (Altshuler et al., 2015; Crandell and Tobalske, 2015), where the wing structure is composed mostly of flight feathers. Hence, understanding the underwater propulsion mechanism of penguins would highlight how they have adapted to aquatic life.
Three-dimensional motion measurements of the body and wings are necessary for quantitative hydrodynamic investigation. In particular, the angle of attack (AoA) of the wing has a large effect on thrust generation. However, only a few attempts have been made to measure the wing and body kinematics of penguins. Clark and Bemis (1979) filmed the forward flapping-wing swimming of seven penguin species in a straight water tank from the side at an aquarium. The 2D wingtip trajectory and the relationship between the wingbeat frequency and swimming speed were investigated. The authors reported that the length-specific speed is positively correlated with the wingbeat frequency. Hui (1988a) conducted a similar 2D motion analysis of the wing by filming forward flapping-wing swimming of Humboldt penguins (Spheniscus humboldti) in a straight water tank from the side with a single cine camera at an aquarium. The results suggested that the thrust is generated by controlling the AoA appropriately during both upstroke and downstroke. Although the 2D motion analyses mentioned above provide a qualitative concept of the lift-based propulsion mechanism (Bannasch, 1995), the 3D relationship between the wing motion, relative flow and resultant hydrodynamic force has not been fully clarified. Another important consideration is the wing bending during flapping, which was reported from previous 2D observations (Clark and Bemis, 1979), because the wing deformation could alter the AoA and resultant thrust.
In this paper, we aim to reveal the wing kinematics and propulsion mechanism in penguins considering wing deformation. For this purpose, we investigated the 3D motion of penguins during forward and horizontal swimming using multiple waterproof video cameras at an aquarium. Gentoo penguins (Pygoscelis papua) were chosen for their relatively high-speed foraging (2.3 m s−1) compared with that of other penguin species (Sato et al., 2007) and long migration, up to 268 km from the colony (Wilson et al., 1998). A series of water tunnel experiments with a 3D printed non-flapping wing was also conducted to obtain its steady wing characteristics, followed by quasi-steady hydrodynamic analysis of the thrust in the swimming penguins. Finally, the propulsion mechanism and effect of wing bending on efficiency were evaluated.
wingbeat amplitude at the wingtip
acceleration of the centre of the body
- aCoB,xb, aCoB,yb, aCoB,zb
acceleration of the centre of the body decomposed into the body coordinate system
coefficient in Bannasch (1995)
angle of attack
acceleration of the centre of the left-wing element j perpendicular to the wing surface
power index in Bannasch (1995)
body length of the penguin
body length of the penguin model in Bannasch (1995)
- BT, PBT
beak tip point and its position
drag coefficient in steady flow
wing element chord length of wing element j
wing element chord length of left-wing element j
lift coefficient in steady flow
mean chord length
centre of the body (the geometric centre of the four body tracking points)
cost of transport
- DF, PDF
dorsal front marker point on the body and its position
direct linear transformation
- DR, PDR
dorsal rear marker point on the body and its position
eccentricity of a spheroid
equation of motion
body added-mass force
- Fbody,am,xb, Fbody,am,yb, Fbody,am,zb
body added-mass force decomposed into the body coordinate system
body drag force
component of the body drag force in the xb direction
fluid force generated by the body
fluid force generated by both wings
fluid force generated by both wings calculated by the equation of motion
fluid force generated by the flat wings calculated by quasi-steady analysis
fluid force generated by the original wings calculated by quasi-steady analysis
component of the fluid force generated by the original wings calculated by quasi-steady analysis in the xb direction
volume force (sum of the gravitational force and buoyancy force)
wing added-mass force
- FLwing,am, FRwing,am
wing added-mass force for each of the left and right wings
wing added-mass force of left-wing element j
wing drag force
- FLwing,drag, FRwing,drag
wing drag force for each of the left and right wings
wing drag force of left-wing element j
wing lift force
- FLwing,lift, FRwing,lift
wing lift force for each of the left and right wings
wing lift force of left-wing element j
wing element number
- k1, k2
- LE, PLE
leading-edge marker point on the wing and its position
- Madd,xb, Madd,yb, Madd,zb
body added mass for each body direction
vector perpendicular to the wing plane of left-wing element j
position of the centre of the body at the start of the wingbeat
body coordinate system at the start of the wingbeat
position of the centre of the body at the end of the wingbeat
provisional body coordinate system
body coordinate system
shoulder coordinate system
world coordinate system
position of the leading-edge marker point after rotation
position of the reference point of wing element j
position of the trailing-edge marker point after rotation
position of the centre of wing element j
- PWL, PWR
hydrodynamic power for each of the left and right wings
position of the wingtip point after rotation
major radius of a spheroid
minor radius of a spheroid
wing element width of wing element j
wing element width of left-wing element j
wing element area of wing element j
wing element area of left-wing element j
- TE, PTE
trailing-edge marker point and its position
- TT, PTT
tail tip point and its position
mean swimming speed (time average of the velocity of the centre of the body during one wingbeat cycle)
volume of the penguin
velocity of the centre of the body
- vCoB,xb, vCoB,yb, vCoB,zb
velocity of the centre of the body decomposed into the body coordinate system
- VF, PVF
ventral front marker point on the body and its position
- VR, PVR
ventral rear marker point on the body and its position
relative flow velocity vector at the centre of left-wing element j
relative flow velocity vector of left-wing element j
relative flow velocity vector at the wingtip
- WB, PWB
wingbase marker point and its position
- WT, PWT
wingtip point and its position
angle of attack
angle of attack of wing element j
angle of attack of left-wing element j
time interval between each frame (1/60 s)
hydrodynamic efficiency of propulsion
kinematic viscosity of sea water
sea water density
MATERIALS AND METHODS
Video recording of swimming penguins and camera calibration
Underwater forward swimming of gentoo penguins [Pygoscelis papua (Forster 1781)] in a water tank at an aquarium was recorded with multiple waterproof video cameras placed inside the water tank; the data acquisition was followed by 3D direct linear transformation (DLT) motion analysis (Abdel-Aziz and Karara, 2015). The recording was conducted at Nagasaki Penguin Aquarium (Nagasaki, Japan) on 7 March 2018 and 4 September 2019 (Fig. 1A,B). Three penguins (individuals 1, 2 and 3) were used as summarized in Table 1. Individual 1 was used in both 2018 and 2019 experiments. For each experiment, individual 1 was treated separately in the analysis, considering the difference in experimental conditions. Thus, the penguin IDs were assigned as 1A, 2, 3 and 1B. Penguins 1A and 2 participated in the 2018 experiment, and penguins 3 and 1B participated in the 2019 experiment. The morphological parameters of each penguin ID are also summarized in Table 1. The body length between the beak tip and the tail feather tip in the swimming form was calculated after the 3D motion analysis was performed (Fig. S1A). All experiments were conducted in accordance with the Animal Experiment Management Regulations at the Tokyo Institute of Technology. The experimental protocol was also approved by the Nagasaki Penguin Aquarium. No penguin was harmed during any part of the study.
The dimensions of the water tank were 14 m in length, 4 m in width and 4 m in water depth. The water temperature was between 17 and 18°C. Underwater penguin swimming was captured with 12 (2018) or 14 (2019) waterproof video cameras (GoPro HERO6 Black and GoPro HERO7 Black, GoPro Inc, USA) at a frame rate of 60 frames s−1. The image resolution was 3840×2160 pixels and the exposure time was 1/480 s. The recording file format was MP4 with HEVC (H.265) codec of lossy compression. We used UHS Speed Class 3 (U3) micro SDXC cards (Evo Plus 64 GB, Samsung Electronics Co., Ltd, Korea) to ensure sufficient writing speed. The cameras were placed inside the water tank as shown in Fig. 1A,B. Since each camera was a stand-alone camera and they were not synchronised with each other, we adjusted the timing by recording flashing lights with all the cameras at the same time in advance. Then, the timing of the recorded videos was adjusted according to the captured flashing. This may cause a slight (at most 1/60 s) timing variation between the cameras, which would lead to an error in 3D reconstruction. However, the error is reduced by using multiple cameras. It was confirmed that the drift of time among the cameras, if any, is less than a frame, even after a whole recording session (which usually lasted 40–60 min).
All the cameras were installed in commercial transparent plastic housings for protection. The cameras on the wall and window were attached to the surface via suction cup mounts. The cameras on the floor were inserted in custom-built 3D-printed plastic casings, each of which was fixed to an aluminium alloy beam with a pre-defined angle. These casings made it easy to remove and re-insert the cameras from them for the replacement of batteries and memory cards. Similar casings were also used for the wall- and window-mounted cameras in the 2019 experiment. The beams were connected to a rectangular aluminium alloy frame on the floor (Fig. 1A,B). Since the calibration chains to which the calibration markers were attached (see below for details) were also connected to this metal frame, the relative positions and angles of the markers for each camera were known. At least half of the cuboid calibration frame was visible from each camera.
Six plastic chains were used to form the cuboid calibration frame for the DLT method. The bottom of each chain was connected to the rectangular aluminium alloy frame on the floor, while the top of each chain was connected to a buoy so that the chain was held vertically (Fig. 1A,B). Four or five orange tapes were attached to each chain. The orange tapes and the buoys are used as reference points of the camera calibration, forming a cuboid calibration frame of 3.6 m in length, 1.5 m in width and 2.0 m (7 March 2018) or 2.5 m (4 September 2019) in height. Commercial motion analysis software (DIPP-motion V/3D, Ditect Co., Ltd., Japan) was used for the camera calibration, with which both the lens distortion and the 3D reconstruction with the DLT method were performed. Although each camera view cannot accommodate all the calibration markers, the software can utilise as many calibration markers as available to reconstruct a large calibration frame (yellow dotted lines, Fig. 1A,B) and complete the camera calibration in a single instance.
The metal framework, plastic chains and buoys remained in the tank during most of the recordings but were removed in some sequences. We confirmed that the posture of the cameras did not change after the frame was removed. The recording duration was 40–60 min, while the target penguin swam freely with other penguins. The thick, yellow chains and white buoys were readily visible, and the penguins generally avoided them. Note that the chains stood still and their motion due to waves or flow was imperceptible, as shown in Movies 1 and 2. Although the chain can move when a penguin directly touches it or swims very close to it, the measurement accuracy was not affected because the images of the chains used for the camera calibration were taken when the chains were still.
Three-dimensional position tracking
Twelve points on a penguin's body were tracked and their 3D coordinates were obtained (Fig. 1C,D). Four markers were attached to the body in the sagittal plane: dorsal front (DF), dorsal rear (DR), ventral front (VF) and ventral rear (VR). Three markers were attached to each wing: wingbase (WB), leading edge (LE) and trailing edge (TE). Left and right wingtips (WTs) were tracked without markers. For the body markers, black or white adhesive tape (Tesa, Beiersdorf AG, Germany) was used (Wilson and Wilson, 1989). For the wing markers, additional strips of orange tape were glued with cyanoacrylate glue (Loctite LBR-005, Henkel, Germany) on the black Tesa tape wrapped around the wing. Since two penguins with markers were recorded at the same time, a piece of orange tape was glued on the top of the DR marker for one of the two penguins for identification. Additional markers were wrapped around the left and right ankles of the penguins in the 2019 experiment, but they were not used in this analysis. All of the tapes were removed after the experiment without leaving any noticeable residue or damage to the plumage. Hereafter, the positions of the points are referred to as Ppoint. For example, PDF represents the position of the DF marker. The body length was obtained as the average distance between the tip of the beak and the tip of the tail during forward straight gliding.
Each tracking point was manually tracked on each captured image and its 3D position was calculated with the same motion analysis software used for the camera calibration. We analysed video sequences while all the tracking points were visible from at least two cameras. In most cases, three or more camera images were available for tracking. Before calculating the kinematic parameters (position, angle, velocity or acceleration), the positional data were processed with a custom MATLAB (2020b and 2021a, MathWorks, Inc., USA) script, where the raw time-series data were smoothed with a 17-point weighted moving average with a generalized Hamming window tuned to be a low-pass filter with a cut-off frequency of 8 Hz and gain of 1 at 0 Hz.
To validate the accuracy of the DLT method, a 1.5 m straight aluminium bar with four tape markers attached every 0.5 m was captured in the same water tank and the positions of the markers and the tips were measured. A diver holding the bar moved back and forth at approximately 2 m s−1 in the calibration frame, and the length between neighbouring markers or tips was calculated and compared with the true value. The calculation was performed every four video frames for a total of 35 frames. The mean error for each length section (0.5 m) was 0.016 m at maximum, which is 3.2% of the section length.
The world coordinate system O-xyz was defined by the calibration frame, as shown in Fig. 1A,B. The long side of the calibration frame was taken as x, the short side as y and the direction opposite to gravity as z. To describe the swimming motion and wing kinematics, we further defined two body-fixed coordinate systems that move with the body, i.e. the body coordinate system Ob-xbybzb fixed to the body centre and the shoulder coordinate system Os-xsyszs fixed to PWB (Fig. 2A). The geometric centre of the four body tracking points (PDF, PDR, PVF and PVR) was chosen as the origin of the body coordinate system Ob. Since the centre of mass of the penguin was unknown, Ob was used as the reference point for the body motion in the world coordinate system for obtaining the swimming trajectory, velocity, or acceleration. Hereafter, we refer to Ob as the centre of the body (CoB). The shoulder coordinate system was created by translating the body coordinate system to the wingbase (PWB) (Fig. 2A). Note that the shoulder coordinate system is fixed to the body, not to the wing. The shoulder coordinate system was used to define the wing kinematics relative to the body.
To define the orientation of the body coordinate system, we first considered what would be an appropriate body-axis. The directions of the vectors based on the ventral body markers and the tail-beak line in glide swimming were slightly different. After some examinations, we concluded the latter to be the more appropriate representation of the body-axis because it is more aligned with the swimming direction during gliding (Fig. 2B). Therefore, we defined the body longitudinal axis as follows. We first selected a short forward gliding sequence for each penguin ID. For the sequence, the tail feather tip (TT, PTT) and the beak tip (BT, PBT) were manually tracked. The resultant vector is the direction of the longitudinal axis. To define the body coordinate system for each time frame in each swimming sequence, we first defined a provisional body coordinate system Ob-xb,pyb,pzb,p based on the body markers (Fig. 2B, top). The direction of the xb,p axis is defined as a vector whose origin is Ob and is parallel to the vector . The direction of the yb,p axis is determined by the cross product of the vectors and . The direction of the zb,p axis is determined by the cross product of the unit vectors for xb,p and yb,p. To obtain the deviation of the xb,p axis from the longitudinal axis, the vector was first projected onto the xb,pzb,p plane (sagittal plane). The angle between this projected longitudinal vector and the xb,p axis was calculated for each frame during forward gliding and its average value was obtained. We call this angle the constant pitch-offset angle ε, whose values were 9.5 deg, 1.9 deg, 5.5 deg and −3.5 deg for penguin IDs 1A, 2, 3 and 1B, respectively. Here, ϵ is positive when the xb,p axis points upward relative to as illustrated in Fig. 2B, top. The body coordinate system Ob-xbybzb was then calculated by rotating the Ob-xb,pyb,pzb,p around the yb,p axis (Fig. 2B, bottom). Note that the plane xbzb represents the sagittal plane and the xb axis represents the longitudinal direction of the body.
Selection of the forward and horizontal swimming cases
The forward and horizontal swimming cases were selected by the moving direction of the CoB for each wingbeat (Fig. 2C). The body coordinate system at the start of the wingbeat was defined as Ob0-xb0yb0zb0 and was fixed to the position of the CoB at the start of the wingbeat, Ob0. The position of the CoB at the end of the wingbeat was defined as Ob1. Subsequently, a swimming direction was determined by the vector . The azimuth angle of the swimming direction is a rotational angle of the vector from the xb0zb0 plane around the zb0 axis measured in the counter-clockwise direction. Hence, a positive azimuth angle indicates a left turn. The inclination angle is an angle between the vector and the plane xb0yb0. Finally, forward swimming was defined as swimming with an absolute azimuth angle of 15 deg or less and an absolute inclination angle of 20 deg or less. We also considered the ascending angle in the world coordinate system, which is an angle between the vector and the horizontal plane. To focus on the basic swimming mechanisms with less complexity associated with buoyancy and gravity, we selected only horizontal swimming where the absolute ascending angles are 20 deg or less.
Analysed wingbeat cases and ensemble averaging
For the analysed penguins, 54 sequences were obtained, where the number of wingbeats in each sequence ranged from 0 to 6. The number of wingbeats was 0 in 16 sequences because either the penguins only glided or the flapping kinematics we could observe was shorter than one complete cycle; therefore, these sequences were excluded. In total, 65 wingbeats from 38 sequences were obtained, where 47 wingbeats involved forward (straight) swimming and 18 are turning. Among the 47 wingbeats of forward swimming, 6 were excluded from the analysis because the absolute ascending angles were larger than 20 deg. In addition, one wingbeat of the penguin ID 2 was excluded because we determined that this wingbeat was in the middle of a manoeuvre from observation of the videos and the flexed reconstructed trajectory. Consequently, 40 wingbeats (8, 7, 17 and 8 wingbeats from the penguin IDs 1A, 2, 3 and 1B) from 30 different sequences were analysed. Each wingbeat began with the start of an upstroke and finished with the end of a downstroke. A downstroke always followed soon after an upstroke without delay (i.e. no gliding while the wings are held upwards). Although both short-term (less than a second) and long-term (more than a second) gliding (i.e. swimming without flapping) were often observed between sporadic wingbeats or sets of continuous wingbeats, the effect of the preceding status for each wingbeat, whether flapping or gliding, was not taken into account in this study. To perform ensemble averaging for the time-series kinematics, the time for each wingbeat was normalized by the respective wingbeat period, followed by a linear interpolation, which generated 50 data points for each wingbeat. Hereafter, we express this normalized time as t′. The mean value and standard deviation (s.d.) were calculated from the ensemble averaging for each penguin ID. Then, the pooled mean value and s.d. of the mean curve across the four penguin IDs were obtained. We show the pooled mean and s.d. in the main text and figures unless otherwise noted. The values are presented as the means±s.d. A P-value less than 0.05 was considered statistically significant.
Calculation of the velocity and acceleration
Wing motions relative to the body
We defined three rotational angles of the wing in the shoulder coordinate system to express the 3D wing kinematics relative to the body as follows: flapping angle, βflap, sweepback angle, βsweep, and feathering angle, βfeather (Fig. 3A). The flapping angle, βflap, which represents the main flapping motion, is an angle between the xsys plane and orthogonal projection of the line segment PWBPLE onto the yszs plane (Fig. 3C). Subsequently, PLE and PTE are rotated by −βflap around the xs axis so that PLE is on the xsys plane. The rotated PLE and PTE are referred to as PLE′ and PTE′, respectively. Then, the sweepback angle, βsweep, is defined as an angle between the line PWBPLE′ and the ys axis; hence, its positive direction indicates that the wing is swept backwards. Note that the βsweep is slightly smaller than the angle between the ys axis and the orthogonal projection of the line PWBPLE onto the xsys plane (Fig. 3B). The feathering angle βfeather is an angle between the rotated inner wing (i.e. a triangle composed of PWB, PLE′, and PTE′) and the xsys plane. When PTE′ is above the xsys plane, as shown in Fig. 3A, βfeather is defined as negative.
We also defined a flapping plane for each wingbeat as a plane that is parallel to the yb axis and contains the maximum and minimum wingtip (PWT) positions during the wingbeat. Our definition of the flapping plane is similar to the ‘anatomical’ stroke plane, which represents the wing trajectory relative to the mid-frontal plane of the bird (Tobalske et al., 2007).
The start and end times of the upstroke and downstroke were determined by the time variation of the flapping angle βflap. Since the wings have anhedral angles (i.e. negative βflap) at a relaxed state (i.e. gliding), as demonstrated in Fig. S1A, the wingbeats always start with an upstroke. Therefore, we chose the beginning of the upstroke as the beginning of the single wingbeat cycle. The start and end times of the upstrokes were defined as the periods when βflap attains its negative and positive peaks, respectively. The start of the succeeding downstroke was identical to the end of the preceding upstroke. The end of the downstroke was determined by the time of the subsequent negative peak of βflap. It was sometimes difficult to determine the start and end of the wingbeat cycle since the negative peaks of βflap were not always obvious when the penguin flapped only once during the recorded sequence. Therefore, when the negative peak of βflap was not obvious, the time when the absolute value of the derivative of βflap exceeded a threshold of 50 deg s−1 was chosen as the start and end of the wingbeat cycle.
Each representative wingbeat period for a pair of wings was determined by averaging the periods of both wings' wingbeats, considering the timing difference between the left and right wings. The wingbeat frequency f was then calculated as 1/(wingbeat period).
Wing deformation and flat wing model
The wing deformation was evaluated by out-of-plane bending and in-plane folding (Fig. 3D). The bending angle θbend is an angle between the inner wing formed by PWB, PLE and PTE and an outer wing formed by PLE, PTE and PWT. A positive θbend indicates upward bending and a negative θbend indicates downward bending. To calculate the in-plane folding, the outer wing was rotated around a line connecting PLE and PTE by −θbend so that the rotated wingtip point (PWT′) lay in the inner-wing plane and the whole wing was flat. Here, we define a flat wing model, which is composed of the original wing's inner wing and the rotated outer wing. Thus, θbend of the flat wing model is zero. The folding angle, θfold, was then calculated as the angle between the inner wing leading edge (the line segment PWBPLE) and the outer wing leading edge of the flat wing model (the line segment PLEPWT′). A positive θfold angle indicates that the wingtip of the flat wing model (PWT′) is located posterior to the line PWBPLE. Note that θfold is generally positive because the leading edge is curved backward towards the wingtip (Fig. 1C).
Equation of motion of swimming penguins
To estimate the volume of the penguin body for calculating the buoyancy force, we created body models for each penguin ID (Fig. S1A). The lateral outline was obtained from a side photograph. The cross-sectional shape was assumed to be circular. To compensate for the lens distortion and refraction at the window of the tank, the aspect ratio of the model was corrected based on the body length and the distance between dorsal and ventral markers.
Quasi-steady calculation of the hydrodynamic force of the body
A 3D wing model was created for morphological measurement and hydrodynamics investigation. The 3D shape of the right wing of individual 1 was measured with a 3D scanner (Space Spider, Artec3D, Luxembourg) on 27 September 2016 (Fig. S1B). Coloured circular stickers were attached to the surface as visual features to improve the scanning quality. In addition, since the wing is a thin object and the scanner easily lost the target when moving between the dorsal and ventral surfaces, a plastic clothes peg with a weak clamping force was attached to the leading edge to provide a suitable silhouette for the scanner. During the scanning, which lasted a few minutes, the penguin was standing on the floor with one of its wings outstretched to the side of the body. The wingtip was gently held by the fingers of the scanner operator while the body was gently held by a staff member of the aquarium. A 3D wing shape was reconstructed from the scanned data with scanning/reconstruction software (Artec Studio 11 Professional, Artec 3D, Luxembourg). The initial data consisted of several sets of 3D point clouds. They were manually integrated using the surface features. The clothes peg was erased at this stage. From the point cloud, a tetrahedron wing surface mesh was constructed with the ‘sharp fusion’ option, which is the most accurate option for mesh generation. To create a watertight wing model with smooth surfaces for the following hydrodynamic experiment, the mesh was exported from Artec Studio and imported to CAD software (Rhinoceros 5, Robert McNeel & Associates, USA), where cross-sectional profile curves were generated and smoothed. The wing torsion (spanwise twist) was small: the maximum and minimum angle between the wing chord and the horizon was 8.4 deg and −7.2 deg, respectively. Note that the torsion could be affected by the holding posture during the scanning. The camber height was small too, ranging between −1.6% and 4.3% (Table S1). Therefore, the cross sections were modified into symmetrical aerofoils having no camber, and the wing chords at different spanwise locations were aligned so there is no spanwise torsion (Fig. S1C,D). This helped to clarify how the wing kinematics and deformation affect the hydrodynamic force in the quasi-steady calculation. The curves were then used to generate a smooth 3D wing surface. The wing length and wing projected area of the created symmetric model were 0.252 m and 0.014 m2, respectively. Note that the scan data and wing model above are the same as those used in our previous work (Shen et al., 2021).
For the quasi-steady calculation of the force on each wing, we decomposed each wing into three elements (Fig. S1E). The outline of the wing model was scaled and aligned to the time-averaged shape of the flat wing model for each penguin so that the distance between each vertex and the outline was minimized. Time averaging was conducted for all the wingbeats to generate a single averaged shape for each wing. The outer wing was bisected by the chord line parallel to the line PLEPTE. The element area, Sj, and element width, rj, were obtained for each element, where subscript j is the wing element number. The chord length, cj, was defined as Sj/rj. The centre of each element, PWC,j, was calculated to be the geometric centre of the vertices of the element. The above operation was repeated for each wing and for each penguin ID.
Quasi-steady calculation of the hydrodynamic force of the wing
By substituting Eqns 21–23 into Eqn 17, we obtain Ffluid,wing,QS, the quasi-steady wing force. Following the same procedure, the force on the flat wing was calculated, denoted by Ffluid,wing,flat,QS.
Fluid dynamics parameters
The Reynolds number (Re) and Strouhal number (St) are fundamental fluid dynamics parameters for the flow past an oscillating object. Re is a dimensionless number that represents the ratio of the inertial effect to the viscous effect of the flow. St is a dimensionless number that represents the unsteadiness of the transverse oscillation of the object relative to the translational main flow. For our penguin with flapping wings, we defined the Reynolds number of the wing as Re=|vWT|cmean/ν, where ν is the kinematic viscosity of seawater at 17°C (i.e. 1.1304×10–6 m2 s−1) (Fresh Water and Seawater Properties: https://ittc.info/media/4048/75-02-01-03.pdf, accessed 7 July 2021), |vWT| is the magnitude of relative flow velocity for the wing at PWT (m s−1) and cmean is the mean chord length, defined as S/r. Here, S is the wing projected area and r is the wing length, i.e. the distance between PWB and PWT. Since the relative flow velocity varies with the wing motion, the value of Re also varies within one wingbeat. The Strouhal number was defined as St=fA/Ub, where f is the wingbeat frequency (Hz), A is the wingbeat amplitude at the wingtip (m) and Ub is the mean swimming speed (m s−1) (similar to Taylor et al., 2003). The wingbeat amplitude was determined from the wingtip trajectory viewed in the xbzb plane (i.e. lateral view of the body) as the length in the zb direction (i.e. the maximum zb minus the minimum zb of the trajectory). The mean swimming speed was defined as the time average of the velocity of the CoB during one wingbeat cycle. The wingbeat frequency was defined as the inverse of a single wingbeat period directly determined by the measured flapping angle. Note that previous biologging studies calculated the wingbeat frequency of penguins as the number of wingbeats per unit time based on the recorded data of the accelerometer for a longer period (Sato et al., 2007, 2010; Watanuki et al., 2006), while our definition focuses on each wingbeat.
Overview of the measured flapping-wing swimming
The measured kinematic parameters for the wingbeats of forward and horizontal swimming are summarized in Table 2. See also Movies 1 and 2 for example sequence and example wingbeat from multiple cropped views, respectively. The ensemble average and s.d. of each parameter in a normalized wingbeat cycle within each penguin ID are summarized in Figs S2 and S3.
The measured mean swimming speed Ub for all four penguin IDs ranged from 0.83 to 2.07 m s−1, and their mean and s.d. values were 1.17±0.13, 1.06±0.42, 1.44±0.30 and 1.45±0.26 m s−1 for the four penguin IDs. The wingbeat frequency ranged from 1.43 to 2.50 Hz, and the mean and s.d. values were 1.84±0.22, 2.01±0.28, 2.09±0.28 and 2.06±0.22 Hz. The mean ratios of downstroke duration to wingbeat period were 0.44±0.07, 0.42±0.10, 0.50±0.06 and 0.49±0.06. Thus, durations of the upstroke and downstroke tended to be similar. For some wingbeats, however, the downstroke durations were notably longer than the upstroke durations (Table S2), suggesting that the penguin is capable of modulating the upstroke and downstroke durations.
Velocity and acceleration of the body
As a result, both upstroke and downstroke produced a forward acceleration in the longitudinal body-axis (aCoB,xb>0). The maximum aCoB,xb in the upstroke (1.21±0.63 m s−2) at t′= 0.26 was larger than that in the downstroke (0.82±0.60 m s−2) at t′=0.76. The time-averaged aCoB,xb during the upstroke of 0.34±0.21 m s−2 was also larger than that during the downstroke (0.20±0.17 m s−2). The maximum net acceleration |aCoB| in the upstroke (2.68±0.71 m s−2) was also larger than that in the downstroke (2.09±0.69 m s−2). Taking the increment of vCoB,xb for the upstroke and downstroke to evaluate the contribution to the gain of the forward velocity, the upstroke resulted in a gain of 0.10±0.07 m s−1 and the downstroke resulted in a gain of 0.06±0.06 m s−1. Thus, both upstroke and downstroke increased the forward velocity, but the upstroke contributed more than the downstroke.
The dorsoventral acceleration aCoB,zb sinusoidally changed in synchronism with the flapping motion, indicating that the body accelerated towards the ventral direction in the upstroke and accelerated in the dorsal direction in the downstroke. The lateral acceleration aCoB,yb was almost zero, validating our classification of forward swimming.
Ensemble-averaged trajectories of the wingtip and wingbase in the body coordinate system projected onto the frontal, lateral and dorsal views are shown in Fig. 5. Although the wingtip trajectory changed every wingbeat, as indicated by the large s.d. (cross hairs in the figure), the mean trajectory in the lateral view demonstrates that the wingtip followed almost the same straight path during the upstroke and downstroke. The flapping plane was tilted forward with the flapping plane angles (an angle between the flapping plane and xb axis in the body coordinate system) of 75.6±7.0 deg, 71.0±9.1 deg, 76.2±6.2 deg and 77.0±6.8 deg for the four penguin IDs. In the dorsal view, the wingtip was always behind the wingbase, indicating that the wing was always swept backwards.
The ensemble-averaged wing angles are shown in Fig. 4C. The flapping angle βflap sinusoidally increased from −37.1±10.3 deg to 35.8±4.9 deg in the upstroke and then decreased back to −39.9±13.2 deg. The mean βflap over the wingbeat cycle was −2.3±6.7 deg. The sweepback angle βsweep also sinusoidally varied, while its phase was opposite to that of the flapping angle. The βsweep decreased from 37.0±3.6 deg to 15.6±1.9 deg in the upstroke and increased to 37.7±4.4 deg in the downstroke. The mean βsweep was 26.0±2.2 deg. The parameter βsweep was positive during most of the wingbeats, indicating that the wing was regarded as a ‘swept wing’ relative to the body. The waveform of the feathering angle βfeather preceded the flapping angle by approximately π/2. The βfeather increased from −6.7±7.4 deg to the positive peak of 6.9±5.9 deg in the upstroke and then decreased from −5.5±2.7 deg to the negative peak of −25.7±4.9 deg in the downstroke. The mean βfeather was −7.1±4.3 deg. Hence, pronation (i.e. negative βfeather) in the downstroke exceeded supination in the upstroke.
The range of Re was between 2.2×104 and 19.0×104 (Table 2B).
Wing deformation and angle of attack
Considerable bending occurred in synchronism with flapping (Fig. 4D). The bending angle θbend started from −4.8±2.7 deg and exhibited a negative peak of −22.7±3.3 deg at t′=0.32 in the upstroke. Then, θbend increased back beyond the initial value, reaching 10.6±3.3 deg at t′=0.80. Therefore, bending was more pronounced in the upstroke than in the downstroke. On the other hand, the folding angle θfold remained almost constant. The mean θfold was 22.0±1.3 deg. Since the original θfold approximated from a top-view photo of the expanded wing (Fig. 1C) was 20 deg, the in-plane deformation can be regarded as negligible.
The AoA varied in synchronism with the flapping, while its magnitude and phase differed with respect to the location and bending (Fig. 4E). The amplitude of the AoA increased and the phase was delayed as the spanwise location of the cross section of interest moved towards the wingtip. The AoAs in the original wings (solid lines in Fig. 4E) were generally negative during the upstroke and positive during the downstroke because of the feathering. The minimum and maximum AoAs were −16.5±3.0 deg and 5.7±3.4 deg for element 1, −17.7±2.0 deg and 17.9±3.8 deg for element 2, and −29.0±3.3 deg and 30.4±1.4 deg for element 3. Comparison with the flat wing model (dashed lines in Fig. 4E) suggests that the wing bending reduced the magnitude of the AoA in the first half of each stroke. In particular, the AoA amplitude at wing element 2 in the upstroke was remarkably reduced by the wing bending. The minimum and maximum AoAs of the flat wing model were −25.6±3.1 deg and 22.1±3.3 deg for element 2 and −32.6±3.3 deg and 33.3±3.0 deg for element 3 (note the element 1 of the flat wing model is the same as that in the original wing model).
Cross-sectional profiles and hydrodynamic characteristics of the wing
The wing profiles measured by 3D scanning were found to be close to the symmetric aerofoils as summarized in Table S1. The ratio of the maximum thickness to the chord length was 18.6% at a chordwise position of 30%. The ratio of the camber height to the chord length was approximately zero. These results are consistent with the previous report for three species of penguins by Bannasch (1995).
The lift, drag and lift-to-drag ratio for the steady case and the cross-flow vortex model are summarized in Fig. 6A,B. The maximum lift coefficient for the steady case CL,steady was 0.83 at α=16.5 deg and the maximum lift-to-drag ratio was 7.7 at α=9.0 deg. The lift slope near AoA=0, (dCL,steady/dα)|limα→0, was 2.80 (rad−1). A stall occurred at approximately α=16.5 deg. After the stall, CL,steady dropped and remained at approximately 0.7 until the AoA reached 40 deg, while the drag coefficient CD,steady abruptly increased just after the stall and continued to increase. In the cross-flow vortex model, the wing does not stall and the basic steady lift and drag characteristics do not change, as explained in the ‘Quasi-steady calculation of the hydrodynamic force of the wing’ in the Materials and Methods. CL and CD of the model matched CL,steady and CD,steady before the stall, respectively.
Quasi-steady calculation of the forces and efficiency
The forces calculated by the quasi-steady method are shown in Fig. 6C–F. In the forward (xb) direction (Fig. 6C), the wing force for the original wing (red) qualitatively matched the value obtained from the EoM (pink), that is, there were positive peaks in both mid-upstroke and mid-downstroke. The drag force of the body (purple), the added-mass force of the body (green) and the volume force (gravity and buoyancy) (light blue) made little contribution to the total force (black) in the xb direction. The virtual suppression of bending (blue) resulted in a loss in the thrust (forward force) during the upstroke. The average thrusts over the upstroke, downstroke and whole cycle for the original wing were 1.86±1.34 N, 5.75±2.28 N and 3.61±1.64 N, while those for the flat wing were 0.80±1.56 N, 5.89±2.56 N and 3.06±1.69 N, respectively. Therefore, wing bending enhanced the thrust during the upstroke, by a factor of 2.3, while the thrust during the downstroke was slightly decreased by a factor of 0.98. The thrust per wingbeat was enhanced by a factor of 1.18 by the wing bending. The hydrodynamic mechanism of propulsion is explained in the Discussion.
The difference between the wing forces calculated by the quasi-steady method and the wing force estimated from the EoM is more prominent when we examine them in terms of the magnitude (Fig. 6F). Note that the magnitudes of the forces were calculated in the same way as in Eqns 29 and 30. Compared with the wing force estimated from the EoM (pink), the quasi-steady wing force (red) was smaller during the upstroke but larger during the downstroke. Notably, the magnitude of the wing force for the flat wing (blue) was much larger than that for the original wing (red), unlike the forward component. This may be because the AoA for the flat wing was larger than that for the original wing, resulting in larger drag due to the smaller lift-to-drag ratio and thus spoiling the thrust (forward-facing) component of the lift.
The hydrodynamic efficiency of propulsion of the original wing during the wingbeat was 0.27±0.06, while that of the flat wing was 0.15±0.07. Therefore, the wing bending can enhance the efficiency by a factor of 1.8. In each penguin ID (summarized in Fig. S4), the efficiency for the original wing was significantly higher than the efficiency for the flat wing, that is, 1.8, 4.2, 1.5 and 1.6 times higher for the four penguin IDs (paired t-test, P<10−5 for the four penguin IDs). The hydrodynamic efficiency of propulsion during the upstroke was 0.12±0.11 for the original wing, while that for the flat wing was 0.02±0.11. That is, the wing bending can enhance the efficiency during the upstroke by a factor of 5.6. On the other hand, the hydrodynamic efficiency of propulsion during the downstroke was 0.54±0.15 for the original wing, while that for the flat wing was 0.37±0.08. Therefore, the improvement of the efficiency during the downstroke was a factor of 1.5, which is less significant than that during upstroke.
Strouhal number with the mean swimming speed
There was a positive correlation between the mean swimming speed and the wingbeat frequency (linear regression, R2=0.173, P=0.008) (Fig. 7A). The wingbeat amplitude was also correlated with the mean swimming speed (linear regression, R2=0.083, P=0.009) (Fig. 7B). Since the frequency and amplitude tended to increase with increasing mean swimming speed, the product of the two, fA, also increased with increasing mean swimming speed (Fig. 7C). However, the increase in the mean swimming speed was more rapid, so the Strouhal number (St=fA/Ub) tended to decrease as the mean swimming speed increased (linear regression, R2=0.294, P<10−8) (Fig. 7D). In particular, when the mean swimming speed was larger than 1.5 m s−1, most of the Strouhal numbers (17 of 20) fell within the range of 0.2 to 0.4, which is expected to be efficient for oscillating-wing propulsion (Anderson et al., 1998; Triantafyllou et al., 1991).
Contribution of the upstroke and downstroke to propulsion
Previous biologging studies reported the velocity and acceleration of various penguins (Sato et al., 2002, 2007; Watanuki et al., 2006), but the details within a wingbeat have never been examined. Our results show that both upstroke and downstroke contributed to the gain of forward velocity (Fig. 4A). This balanced gain of velocity is expected to lead to efficient swimming (Lovvorn and Liggins, 2002). As for the force magnitude, both the quasi-steady calculation Ffluid,wing,QS and equation of motion Ffluid,wing,EoM showed that the upstroke generated larger force than the downstroke (Fig. 6F). This is probably because dive depth in our measurement was shallow: approximately 1–3 m. Swimming at a shallow depth requires net downward force to balance the upward buoyancy. This in turn favours greater force production during the upstroke than the downstroke. At a deep depth with smaller upward buoyancy due to compression of air volume, the wing kinematics may deviate from that at a shallow depth as in our measurements.
A previous study of wild penguins reported that both upstroke and downstroke caused a surging acceleration: Watanuki et al. (2006) measured the surge (anterior–posterior) and heave (dorsoventral) accelerations of little penguins (Eudyptula minor) while descending swimming using accelerometers. The maximum surge acceleration during the upstroke was 2.36 m s−2, which was 20% greater than the acceleration of 1.97 m s−2 during the downstroke. The maximum magnitude of the heave acceleration during the downstroke was 4.94 m s−2, which was 10% greater than that of 4.48 m s−2 during the upstroke. Hence, our study supports that the contribution of both upstroke and downstroke to propulsion is common in penguins.
Stroke-acceleration patterns have also attracted attention in other wing-propelled diving animals such as alcids, sea turtles and sea lions. Previous studies of alcids revealed that alcids also accelerate forward during both upstroke and downstroke, while the downstroke produces notably larger acceleration than the upstroke (Johansson and Aldrin, 2002; Lapsansky and Tobalske, 2019; Watanuki et al., 2006) unlike our penguins. Previous kinematic studies of sea turtles reported that upstroke generates little acceleration and downstroke generates more acceleration (Becking et al., 2004; Booth, 2014; Davenport et al., 1984; Rivera et al., 2009) unlike our penguins or the alcids. For sea lions, although no direct measurements of acceleration have been made, the results of a kinematic study suggested that a large amount of thrust is generated at the end of the downstroke (‘paddle phase’; Feldkamp, 1987). In summary, penguins are the only wing-propelled diving animal whose upstroke contributes to propulsion as much as, or more than, the downstroke. Penguins have a pair of particularly large pectoral muscles for upstroke (supracoracoideuses) compared with other birds (Baldwin, 1988; Bannasch, 1994), to which the strong upstroke is attributed.
Thrust generation and wing deformation
Our 3D measurement of the wing and body kinematics makes it possible to investigate the hydrodynamic mechanism of thrust generation and the effect of wing deformation via a quasi-steady calculation. Similarly to a previous description by Azuma (2006), the flapping plane of our penguin is tilted forward; i.e. the wingtip at the beginning of the downstroke is forward (anterior) of the wingtip at the end of the downstroke (Figs 5 and 8), which is opposite to the case of many flying birds whose flapping planes (stroke planes) are tilted backwards in cruising flight (e.g. Tobalske and Dial, 1996; Tobalske et al., 2007). The forward tilting of the flapping plane in penguins has been assumed to produce a downward force that balances the upward buoyant force by lift-based force generation (Clark and Bemis, 1979; Hui, 1988a), as illustrated in Fig. 8A–C. Our quasi-steady calculation confirms this mechanism, as shown in Fig. 8D–F (also see Figs 6 and S3). The upstroke generates wing force in the obliquely downward direction with respect to the normal direction of the flapping plane, which is attributed to the negative AoA modulated by feathering (Fig. 4C,E). Similarly, the downstroke generates wing force in the obliquely upward direction owing to the positive AoA. The resultant mean wing force for the entire wingbeat cycle points to the obliquely downward direction with respect to the normal direction of the flapping plane because the wing force magnitude is greater in upstroke than in downstroke.
Moreover, the calculated wing force provides insight into the source of wing deformation. The peaks of the bending angle are found at the mid-upstroke and mid-downstroke (Fig. 4D), where the hydrodynamic force magnitude also marks the peak values (Fig. 6F). On the other hand, the inertia and added-mass forces are the smallest at the mid-strokes, as the acceleration is small. Therefore, it is reasonable to assume that the bending deformation of the wing is passively caused by the hydrodynamic force rather than the inertial force.
Note that the quasi-steady analysis is sensitive to the lift and drag coefficients and the AoA used in the calculation. The actual CL and CD could differ from the present CL and CD generated by the cross-flow vortex model assuming that the wing does not stall during flapping (Bandyopadhyay et al., 2008). Nevertheless, our quasi-steady calculation can visualize the relationship between the wing kinematics, AoA and resultant wing force. Future studies may explore the impact of the shape of the CL and CD curves on the estimated forces.
Torsional deformation of the wing could not be quantified in this study due to the limitation in the number of wing markers. A recent anatomical study reported that the trailing edge region of a penguin flipper is relatively flexible because the skeleton does not extend to the region (DeBlois and Motani, 2019), implying the possibility of chordwise bending during flapping. Without torsion or chordwise bending, the AoA during flapping increases towards the wingtip where the moving velocity is greatest. Therefore, it is also possible that penguins properly utilise torsional or chordwise deformation to suppress excessive AoA near the wingtip to generate more thrust, as in ‘washout’ in flying birds (Deetjen et al., 2017; Maeda et al., 2017; Wolf and Konrath, 2015) or insects (e.g. Walker et al., 2009, 2010; Zheng et al., 2013). Inclusion of torsional deformation may explain some of the difference in wing force between Ffluid,wing,QS and Ffluid,wing,EoM, particularly during the downstroke (Fig. 6C–F). Further biological measurements and hydrodynamic investigations are desired to clarify the advantages and disadvantages of deformable wings compared with rigid wings.
Dorsoventral and lateral forces
The lift-based flapping propulsion is accompanied by dorsoventral and lateral forces. Inefficiency arises if these dorsoventral and lateral forces are large compared with the thrust. The quasi-steady calculation revealed that the dorsoventral (zb) force is significant (Fig. 8D,E; see also Fig. 6E and Fig. S3D,I,N,S). Although the lateral (yb) forces by a pair of wings are cancelled out (Fig. 6D and Fig. S3B,G,L,Q), the lateral force by each wing is substantial (Fig. S3C,H,M,R). Notably, elimination of the wing bending increases the wing forces in both zb and single-wing yb. The pooled mean of the absolute values in zb for the flat wings (25.16±5.06 N) is larger than that of the original wings (19.32±4.28 N). The pooled mean of the absolute values in yb for the single flat wing (6.05±1.86 N) is larger than that of the single original wing (4.52±1.59 N). Thus, wing bending may save energy by reducing the dorsoventral and lateral forces.
The large dorsoventral force induces oscillation of the body (vCoB,zb and aCoB,zb in Fig. 4A,B). The body oscillation perpendicular to the main flow may increase both friction drag due to thinning of the boundary layer (Lighthill, 1971) and pressure drag due to an increase in the AoA of the body. The positive and negative peak values of vCoB,zb were 0.14±0.07 m s−1 and −0.12±0.04 m s−1, respectively (Fig. 4A). The minimum, maximum and mean values of the AoA of the body (the angle between xb axis and the orthogonal projection of the vCoB onto the xbzb plane) were −7.7±2.2 deg, 6.7±1.6 deg and −1.0±1.7 deg, respectively.
For manoeuvres such as turning, ascending or descending, these large dorsoventral and lateral wing forces may be beneficial if the AoA of the wing can be changed for each half-stroke. For example, if the AoA during a downstroke is set to be zero without changing the preceding upstroke kinematics, the penguin would translate to its ventral direction owing to the downward force produced in the preceding upstroke. Similarly, if the AoA of the left wing were set to zero without changing the right wing kinematics, the penguin would translate to the left due mainly to the lateral force of the right wing (and the body would yaw left owing to the differences in forward force). To further investigate the manoeuvrability, motion measurement and hydrodynamic analysis of manoeuvring swimming are needed.
Efficiency of propulsion and the effect of wing deformation
The lift-based thrust generation with oscillating wings (or fins) is employed by various aquatic vertebrate species, such as high-speed fish, cetaceans, sea lions and seals, as summarized by Fish (1996). The efficiency of propulsion has been analysed hydrodynamically and metabolically for various animals. One of the important factors associated with the efficiency of propulsion is the Strouhal number (St). The idea of using this number with swimming organisms originated from Triantafyllou et al. (1991, 1993). The St value represents how often vortices are created and how close they are (Triantafyllou and Triantafyllou, 1995). It is well known that the St values for the oscillating wings of flying or swimming animals fall within a certain range (Taylor et al., 2003). Hydrodynamic experiments on heaving or flapping wings suggest that propulsive efficiency (the ratio of thrust power to input power) is maximized when St is between 0.2 and 0.4 (Anderson et al., 1998; Triantafyllou et al., 1991). For example, the relationship between St and propulsive efficiency has been investigated for cetaceans (Rohr and Fish, 2004) and mantas (Fish et al., 2016, 2017).
As far as we know, we are the first to report St for penguins, which was made possible by our measurement of the wingbeat amplitude (see Materials and Methods for our definition of St). St tends to decrease as the mean swimming speed increases (Fig. 7D) because the change in the wingbeat frequency and amplitude with the speed is small (Fig. 7A,B). Taylor et al. (2003) suggested that the wingbeat frequency and amplitude are presumably restricted by physiology and morphology, thus the St will strongly depend on the speed. Our St falls into the range of 0.2 to 0.4 when the speed is greater than 1.5 m s−1 (Fig. 7D), suggesting that the wing kinematics may be optimised for relatively high-speed swimming of more than 1.5 m s−1.
The swimming efficiency of penguins was previously estimated by the cost of transport (CoT), which is the energy required to move a unit mass for a unit distance. Data from biological measurements of the metabolic rate and the travel distance were used in the calculations. Swimming speed and energy expenditure of various species of penguins during resting, surface or underwater swimming were measured in the swim canal with a respiration chamber, followed by CoT calculation for a wide range of swimming speeds (Baudinette and Gill, 1985; Culik and Wilson, 1991; Culik et al., 1991, 1994; Hui, 1988b). In these studies, it was suggested that the optimal swimming speed at which the CoT is minimized is approximately 2 m s−1. Based on our measurements, St at 2 m s−1 is expected to be in the range of 0.2 to 0.4 (Fig. 7D). In addition, the foraging speeds of various species of penguins in the wild measured with a propeller-type velocimeter were also approximately 2 m s−1 (Sato et al., 2007). These studies are in good agreement in terms of the efficiency of penguin swimming.
In this paper, we defined the efficiency as in Eqn 28, and the original wing model and flat wing model were compared. As a result, we found that the original wing has higher efficiency than the flat wing (Fig. S4). This enhancement of the efficiency is assumed to be due to the reduction in dorsoventral (zb) and lateral (yb) forces of each wing by the wing bending. Therefore, the bending deformation of the wings may reduce the non-propulsive force and improve propulsive efficiency in penguins. As for the comparison with other animal taxa propelling with flapping wings, Fish et al. (2016) examined the swimming strokes of the manta (Manta birostris) based on kinematic analysis and computational fluid dynamic models. They claimed that the bending deformation of the flapping fin is associated with high propulsive efficiency, perhaps owing to reduction in induced drag. On the other hand, our quasi-steady calculation cannot capture the details of unsteady flow structure around the wing and resultant change in induced drag. To investigate this point, hydrodynamic experiment or simulation of unsteady flapping wing is desired.
Anatomical movability of the shoulder, elbow, wrist and digit joints
Several anatomical studies have been conducted on the wing structure of penguins. Raikow et al. (1988) measured the range of flexion (‘folding’ in our study) at the elbow, wrist and major digit joints using the carcasses of penguins, alcids, diving-petrels and non-diving birds. They showed that the range of folding angle at those joints in penguins were substantially lower than those for other species. Since the LE marker in our experiment is located near the wrist joint, our folding angle (θfold) is an indicator of how much folding occurs at the wrist and major digit joints. Our θfold remained at approximately 20 deg during a wingbeat (Fig. 4D), while the sum of the flexion range of the wrist and major digit joints is 54 deg according to Raikow et al. (1988). These results suggest that the penguin's wrist and hand are fully extended and almost fixed during wingbeats even though they can flex in structure. On the other hand, neither the bending range nor the bending stiffness of the penguin wing has been reported so far. Our measurements showed that the bending deformation occurs during the wingbeat (Fig. 4D). Moreover, this bending affects the AoA of the wing (Fig. 4E).
Lucas et al. (2014) suggested that wing and fin bending during flight or swimming in various species generally occurs at 60–70% of wing length from the wingbase. By contrast, the wing bending in our model is calculated at the PLE–PTE line located at 30–40% of the wing length from the wingbase. Although penguins are not included in the work of Lucas et al., the same may be true for penguins because of the possible movability at the major digits. If the bending mainly occurred at the major digits instead of the wrist, our calculation of the bending angle would underestimate the local bending at the major digits. Note that the bending angle (‘flexion angle’ in Lucas et al., 2014) ranged from 20 deg to 40 deg across many species, while the maximum bending was approximately 20 deg in our calculation.
Feathering rotation may occur at the elbow, wrist, or major digit joints; that is, torsional deformation may exist. The torsional deformation is expected to affect the AoA of the wing and resultant wing forces. We therefore need to study the anatomical movability of the wing joints in further detail.
Bannasch (1994) investigated the movability at the shoulder joint of the penguin using carcasses. It was reported that the maximum pronation (negative βfeather in our study) and supination (positive βfeather) were −53 deg and 39 deg, respectively. The maximum elevation (positive βflap in our study) was 45 deg, whereas no restriction was found for depression (negative βflap). The range of wing motion at the shoulder during wingbeats in slow swimming in our measurements (βfeather: −25.7 deg to 6.9 deg; βflap: −39.9 deg to 36.2 deg) (Fig. 4C) does not seem to reach the above structural limitation, suggesting that there is still sufficient room for various manoeuvres.
Effect of the preceding status of the wingbeat
The details of the time-varying data differ with each wingbeat. In particular, the preceding status (flapping or gliding) seems to affect the kinematic and hydrodynamic patterns of each wingbeat. For example, the waveform of the AoA in penguin ID 2 was different from that of the other penguins at the start of the upstroke (Fig. S2J). This is assumed to be because most of the wingbeats of penguin ID 2 involved continuous flapping: the preceding status of 6 out of 7 wingbeats was ‘flapping’ and that of only 1 wingbeat was ‘gliding’. In the case of the continuous wingbeat, the phase of the AoA near the wingtip is likely to delay because of the preceding wing bending. To clarify the difference between the continuous and single-shot wingbeats, further measurements and analysis of various swimming sequences are desired.
Conclusions and future studies
Our study provides basic information on the kinematics of living penguins during flapping swimming. The 3D positions at various locations on the body and wing were obtained with high accuracy compared with those in previous studies (Clark and Bemis, 1979; Hui, 1988a). Moreover, the folding and bending deformation of the wing during swimming was evaluated for the first time. The kinematic analysis revealed that both upstroke and downstroke contribute to increasing the forward speed. The bending deformation of the wing during upstroke was larger than that during the downstroke, while the folding deformation was not obvious. The bending is large enough to affect the AoA of the wing. By using the lift and drag coefficients obtained from water tunnel experiments, the hydrodynamic force generated by the measured wingbeat was estimated by a quasi-steady calculation. The hydrodynamic efficiency of propulsion was also estimated, showing that wing bending improves efficiency.
In future studies, the torsional deformation of the wing needs to be evaluated. The unsteady hydrodynamics of penguins' flapping wings also needs to be investigated more rigorously. Other types of swimming, such as high-speed propulsion and agile manoeuvring should be measured. Variation in the hydrodynamic mechanism among other penguin species is still unknown, so the present kinematic data for forward slow swimming would serve as a basis for future studies on more complex subjects.
We sincerely thank the staff of the Nagasaki Penguin Aquarium for their warm and enthusiastic support for the video recording and 3D scanning of the penguins. We also thank our laboratory members Yusuke Iwasaki for contributing to the early development of the study, including the preliminary data collection and analysis; Akio Kawahara for the calibration frame rehearsal in the Tokyo Institute of Technology swimming pool; Rei Kitahara, Yoshiya Hayakawa, Mizuki Yoshida and Genki Taguchi for helping with the video recording and analysis; Hiroki Kayasuga for helping in the early development of the recording experiment and water tunnel setup; Tatsushi Yamada for creating the wing model; and Sho Katagiri for the force measurements in the water tunnel experiments. We also thank the anonymous reviewers for improving the manuscript.
Conceptualization: H.T., N.H.; Formal analysis: N.H., T.O.; Funding acquisition: H.T.; Investigation: N.H., T.O., M.M., Y.S., H.T.; Methodology: N.H., T.O., M.M., H.T., Y.S., D.M.K.; Supervision: H.T.; Visualization: H.N.; Writing – original draft: H.N., H.T., M.M.; Writing – review & editing: N.H., H.T., M.M., D.M.K., T.O., Y.S.
This work was supported by a KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas ‘Science of Soft Robot’ project funded by the Japan Society for the Promotion of Science under Grant Number JP18H05468.
The authors declare no competing or financial interests.