Kinematics and hydrodynamics analyses of swimming penguins: wing bending improves propulsion performance

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⁻¹ (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⁻¹) 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.

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⁻¹. 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.

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 P_point. For example, P_DF 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⁻¹ 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 O_b-x_by_bz_b fixed to the body centre and the shoulder coordinate system O_s-x_sy_sz_s fixed to P_WB (Fig. 2A). The geometric centre of the four body tracking points (P_DF, P_DR, P_VF and P_VR) was chosen as the origin of the body coordinate system O_b. Since the centre of mass of the penguin was unknown, O_b 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 O_b as the centre of the body (CoB). The shoulder coordinate system was created by translating the body coordinate system to the wingbase (P_WB) (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, P_TT) and the beak tip (BT, P_BT) 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 O_b-x_b,py_b,pz_b,p based on the body markers (Fig. 2B, top). The direction of the x_b,p axis is defined as a vector whose origin is O_b and is parallel to the vector ⁠. The direction of the y_b,p axis is determined by the cross product of the vectors and ⁠. The direction of the z_b,p axis is determined by the cross product of the unit vectors for x_b,p and y_b,p. To obtain the deviation of the x_b,p axis from the longitudinal axis, the vector was first projected onto the x_b,pz_b,p plane (sagittal plane). The angle between this projected longitudinal vector and the x_b,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 x_b,p axis points upward relative to as illustrated in Fig. 2B, top. The body coordinate system O_b-x_by_bz_b was then calculated by rotating the O_b-x_b,py_b,pz_b,p around the y_b,p axis (Fig. 2B, bottom). Note that the plane x_bz_b represents the sagittal plane and the x_b axis represents the longitudinal direction of the body.

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 O_b0-x_b0y_b0z_b0 and was fixed to the position of the CoB at the start of the wingbeat, O_b0. The position of the CoB at the end of the wingbeat was defined as O_b1. 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 x_b0z_b0 plane around the z_b0 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 x_b0y_b0. 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.

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.

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 x_sy_s plane and orthogonal projection of the line segment P_WBP_LE onto the y_sz_s plane (Fig. 3C). Subsequently, P_LE and P_TE are rotated by −β_flap around the x_s axis so that P_LE is on the x_sy_s plane. The rotated P_LE and P_TE are referred to as P_LE′ and P_TE′, respectively. Then, the sweepback angle, β_sweep, is defined as an angle between the line P_WBP_LE′ and the y_s axis; hence, its positive direction indicates that the wing is swept backwards. Note that the β_sweep is slightly smaller than the angle between the y_s axis and the orthogonal projection of the line P_WBP_LE onto the x_sy_s plane (Fig. 3B). The feathering angle β_feather is an angle between the rotated inner wing (i.e. a triangle composed of P_WB, P_LE′, and P_TE′) and the x_sy_s plane. When P_TE′ is above the x_sy_s 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 y_b axis and contains the maximum and minimum wingtip (P_WT) 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⁻¹ 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).

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 P_WB, P_LE and P_TE and an outer wing formed by P_LE, P_TE and P_WT. 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 P_LE and P_TE by −θ_bend so that the rotated wingtip point (P_WT′) 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 P_WBP_LE) and the outer wing leading edge of the flat wing model (the line segment P_LEP_WT′). A positive θ_fold angle indicates that the wingtip of the flat wing model (P_WT′) is located posterior to the line P_WBP_LE. Note that θ_fold is generally positive because the leading edge is curved backward towards the wingtip (Fig. 1C).

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.

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 m², 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 P_LEP_TE. The element area, S_j, and element width, r_j, were obtained for each element, where subscript j is the wing element number. The chord length, c_j, was defined as S_j/r_j. The centre of each element, P_WC,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.

By substituting Eqns 21–23 into Eqn 17, we obtain F_{fluid,wing,QS}, the quasi-steady wing force. Following the same procedure, the force on the flat wing was calculated, denoted by F_{fluid,wing,flat,QS}.

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=|v_WT|c_mean/ν, where ν is the kinematic viscosity of seawater at 17°C (i.e. 1.1304×10^–6 m² s⁻¹) (Fresh Water and Seawater Properties: https://ittc.info/media/4048/75-02-01-03.pdf, accessed 7 July 2021), |v_WT| is the magnitude of relative flow velocity for the wing at P_WT (m s⁻¹) and c_mean 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 P_WB and P_WT. 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/U_b, where f is the wingbeat frequency (Hz), A is the wingbeat amplitude at the wingtip (m) and U_b is the mean swimming speed (m s⁻¹) (similar to Taylor et al., 2003). The wingbeat amplitude was determined from the wingtip trajectory viewed in the x_bz_b plane (i.e. lateral view of the body) as the length in the z_b direction (i.e. the maximum z_b minus the minimum z_b 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.

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 U_b for all four penguin IDs ranged from 0.83 to 2.07 m s⁻¹, 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⁻¹ 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.

As a result, both upstroke and downstroke produced a forward acceleration in the longitudinal body-axis (a_CoB,xb>0). The maximum a_CoB,xb in the upstroke (1.21±0.63 m s⁻²) at t′= 0.26 was larger than that in the downstroke (0.82±0.60 m s⁻²) at t′=0.76. The time-averaged a_CoB,xb during the upstroke of 0.34±0.21 m s⁻² was also larger than that during the downstroke (0.20±0.17 m s⁻²). The maximum net acceleration |a_CoB| in the upstroke (2.68±0.71 m s⁻²) was also larger than that in the downstroke (2.09±0.69 m s⁻²). Taking the increment of v_CoB,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⁻¹ and the downstroke resulted in a gain of 0.06±0.06 m s⁻¹. Thus, both upstroke and downstroke increased the forward velocity, but the upstroke contributed more than the downstroke.

The dorsoventral acceleration a_CoB,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 a_CoB,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 x_b 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×10⁴ and 19.0×10⁴ (Table 2B).

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).

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 C_L,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, (dC_L,steady/dα)|_limα→0, was 2.80 (rad⁻¹). A stall occurred at approximately α=16.5 deg. After the stall, C_L,steady dropped and remained at approximately 0.7 until the AoA reached 40 deg, while the drag coefficient C_D,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. C_L and C_D of the model matched C_L,steady and C_D,steady before the stall, respectively.

The forces calculated by the quasi-steady method are shown in Fig. 6C–F. In the forward (x_b) 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 x_b 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 forces were small in the lateral (y_b) direction (Fig. 6D). In contrast, the wing forces in the dorsoventral (z_b) direction (Fig. 6E) largely changed in response to each stroke.

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⁻⁵ 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.

There was a positive correlation between the mean swimming speed and the wingbeat frequency (linear regression, R²=0.173, P=0.008) (Fig. 7A). The wingbeat amplitude was also correlated with the mean swimming speed (linear regression, R²=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/U_b) tended to decrease as the mean swimming speed increased (linear regression, R²=0.294, P<10⁻⁸) (Fig. 7D). In particular, when the mean swimming speed was larger than 1.5 m s⁻¹, 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).

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 F_{fluid,wing,QS} and equation of motion F_{fluid,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⁻², which was 20% greater than the acceleration of 1.97 m s⁻² during the downstroke. The maximum magnitude of the heave acceleration during the downstroke was 4.94 m s⁻², which was 10% greater than that of 4.48 m s⁻² 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.

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 C_L and C_D could differ from the present C_L and C_D 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 C_L and C_D 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 F_{fluid,wing,QS} and F_{fluid,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.

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 (z_b) force is significant (Fig. 8D,E; see also Fig. 6E and Fig. S3D,I,N,S). Although the lateral (y_b) 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 z_b and single-wing y_b. The pooled mean of the absolute values in z_b 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 y_b 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 (v_CoB,zb and a_CoB,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 v_CoB,zb were 0.14±0.07 m s⁻¹ and −0.12±0.04 m s⁻¹, respectively (Fig. 4A). The minimum, maximum and mean values of the AoA of the body (the angle between x_b axis and the orthogonal projection of the v_CoB onto the x_bz_b 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.

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⁻¹ (Fig. 7D), suggesting that the wing kinematics may be optimised for relatively high-speed swimming of more than 1.5 m s⁻¹.

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⁻¹. Based on our measurements, St at 2 m s⁻¹ 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⁻¹ (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 (z_b) and lateral (y_b) 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.

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 P_LE–P_TE 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.

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.

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.