Honeybees modify flight trajectories in turbulent wind

In windy conditions, the air is turbulent. The strong and intermittent velocity variations of turbulence are invisible to flying animals. Nevertheless, flying animals, not much larger than the smallest scales of turbulence, manage to maneuver these highly fluctuating conditions quite well. Here we quantify honeybee flight with time-resolved three-dimensional tracking in calm conditions and controlled turbulent winds. We find that honeybee mean speed and acceleration are only weakly correlated with the strength of turbulence. In flight, honeybees accelerate slowly and decelerate rapidly, i.e., they break suddenly during turns and then accelerate again. While this behavior is observed in both calm and turbulent conditions, it is increasingly dominant under turbulent conditions where short straight trajectories are broken by turns and increased maneuvering. This flight-crash behavior is reminiscent of turbulence itself. Our observations may help the development of flight strategies for miniature flying robotics under turbulent conditions.


Introduction
How animals manage to navigate complex turbulent flows is important and of interest to biologists, physicists, engineers, and also the general public. For example, mosquito's prefer not to fly in turbulent conditions [1], possibly associated with the high accelerations in the flow [2,3]. Recent studies have investigated the behavior of small organisms such as planktonic gastropod larvae in turbulence [4] or of larger birds such as a golden eagle in atmospheric turbulence [5]. The ability for animals and especially insects to navigate turbulent flows while in flight has been studied to determine how they are capable of maintaining flight stability and control [6]. Such studies can further improve small-scale robotics for visually guided flight [7], energy harvesting strategies [8] and enhance flapping wing designs [9].
The study of insect flight has gathered much interest due to their relative ease of handling and performing controlled laboratory experiments in different environments, such as the study of fruit flies [10,11] and midges [12]. Particular insects of interest have been the honeybee, bumblebee and orchid bee since they are vital to the environment as natural pollinators [13,14]. The majority of experiments studying bee flight are conducted in laboratory settings with field experiments mostly focusing on bee hive entrance activity [15,16], and monitoring for counting statistics [17][18][19][20][21]. Most experiments performed in the laboratory examine bee flight (bumblebees, orchid bees, etc) in laminar flow and study wing aerodynamics and motion [22]. Other experiments have investigated how bee flight and behavior changes due to external factors such as the presence of predators [23]. Experimental studies examining the influence of turbulent flow in laboratory settings investigated changes in flight in the form of how wing dynamics and bee orientations changed while flying in different conditions [24][25][26]. Studies have shown that bees have increased rolling instabilities in turbulent conditions [24,26], with field measurements also showing stability enhancing behavior to improve rolling stability in the form of bees extending their hind legs during flight [27]. Numerical simulations have shown that there are no significant changes in cycle-averaged aerodynamic forces, moments, and flight power in strong turbulent flow conditions as compared to flight in laminar flow. However, the variance of aerodynamic measures increases with turbulence intensity [28].
The natural flight environment of honeybees are turbulent winds with structures spanning a wide range of scales. Atmospheric turbulence consists of time scale of less than 1 s to typically 1 h and the corresponding length scales are from millimeter up to the boundary layer thickness of hundreds of meters [29,30]. In this study we report data on the flight of honeybees in their natural habitat while flying towards and away from the hive. We reconstruct honeybee trajectories and study honeybee flight paths in natural windless conditions and in the presence of artificially generated turbulent wind using a fan and active grid system.

Methods
Experiments were performed at the Max Planck institute for Dynamics and Self-Organization (MPI-DS) in Göttingen, Germany. Experiments were carried out between 7-13 September 2021, between the hours of 11:00-15:00 each day. The weather during the experiments was calm with little to no natural wind and no precipitation during experiments. As shown in figure 1 experiments were performed near the entrance of a hive with the honeybee species Apis mellifera carnica. Our institute at the time of experiments was in possession of three hives and experiments were performed near the entrance of one of these hives.

Turbulence generation and characterization
Wind and turbulence was produced using a fan array and active grid system identical to the one used in wind-tunnel experiments at the variable density wind-tunnel in Göttingen [31,32]. The fan array consisted of 4 fans with three adjustable speed settings that blew air into to the active grid. The active grid was controlled with an in-house program capable of moving each grid flap independently allowing for adjustable turbulence length and time scales, more details on the controlling code and active grid can be found in [33]. Values for few of the flow quantities measured are given in table 1 along with the number of tracks in those conditions.
To identify the characteristics of the turbulent flow we used a thermal anemometer system previously used in high Reynolds number flows [32]. We also used a sonic anemometer used in atmospheric measurements [34] to measure 3D wind speeds. Here we estimate the turbulence intensity and characteristics using data from the thermal anemometer.
The sonic anemometer was used to determine wind direction and speed for when fans were active and for when fans were off [35]. The sonic anemometer had a sampling frequency of 30 Hz. From the data we clearly see the times in which the fans were turned on and off. When fans were on, the wind direction was strongly in the same direction as where the fans were blowing. When fans were off, there was no clear Table 1. Turbulence characteristics for experiments along with the number of tracks in each experiment. Here R λ = u λ/ν is the Taylor microscale Reynolds number, where u is the root mean square fluctuating velocity, λ = 15ν/εu is the Taylor microscale, ν is the kinematic viscosity of the fluid (1.5 × 10 5 m 2 s −1 for air), and ε is the turbulent energy dissipation rate. The mean flow velocity is given by U. The Kolmogorov timescale is τ = (ν/ε) 1/2 and the Kolmogorov lengthscale is η = (ν 3 /ε) 1/4 . Experiment 0 is for when fans were off and experiments 1-8 are for when fans were on and wind was generated. A higher number of tracks were recorded for the first day of experiments on 09/07/2021, Exp. 1 and 2, due to movement and experiment setup around the hive which caused agitation and higher than usual activity at the hive. For the remaining days the honeybees were less agitated and had adjusted to the changes and activity near the hive. The mean track length across all experiments was approximately 40 frames with the longest track lasting 282 frames. direction to the wind and the mean wind speed was 15% of the mean wind speed when the fans were on. This shows that the weather during the experiments was mild with no strong natural winds present. We also checked wind speeds recorded by a cup anemometer on a nearby building at our institute. The cup anemometer also recorded low wind speeds on the days of experiments with slightly higher values recorded than the sonic anemometer. This is most likely due to the fact that the cup anemometer was on the roof of a building and recorded slightly higher wind speeds because of it.

Imaging and tracking
Imaging was done using three GoPro Hero 9 Black cameras recording at 120 frames per second at a resolution of 2704 × 1520. Using GoPro cameras is a relatively simple and cost effective method of performing biological experiments in the field [36]. The cameras were on the linear setting which digitally corrects for optical distortions caused by the lens with internal software. Calibrations were performed using a 51 × 31 grid of circular black dots on a white background that almost spanned the complete imaging volume. The dots were 0.5cm in diameter and 1cm apart from each other. A blank white background was used for shadow imaging of the honeybees similar to methods previously used capable of highly precise imaging and tracking of particles in turbulence [37,38]. The three cameras had different viewing angles with an overlap region which created a measurement volume of approximately 80 × 60 × 60 cm in size. The existing fencing and an additional screen, not within view of the cameras, encourage the honeybees to fly through the measurement volume (see figure 1). The cameras are triggered by a wireless remote via Bluetooth/WiFi that is able to control all three cameras simultaneously. Each experiment in table 1 consists of 6, 5 minute-long videos, totaling to 30 min of video recordings for each of the 9 experiments with different flow characteristics. The limiting factor for imaging with the GoPro cameras used here is firstly the battery life which allows for approximately 1 h of recording. After the batteries are drained, a replacement set is used and a new set of calibration images are acquired. The second limitation is the overheating of the cameras which occurs due to continuous recording. During wind experiments this issue is resolved due to the cooling nature of the blowing wind. However during no wind experiments cameras switch off after overheating and a waiting period is required before resuming experiments. After the waiting period is done and cameras are turned back on, a new set of calibration images are also acquired to account for any changes in camera position during the pressing of the on button located on the camera bodies.
Identifying honeybees from the images and obtaining tracks was done with in-house codes based on previously used methods [37,38]. The raw images were cropped to only include sections with the white background. The images are then background subtracted, inverted, and binarized with an appropriate threshold. The honeybees now appear as bright pixels with a dark background where the honeybee positions are then identified using the center of mass of the bright pixels. The positions of the honeybees from the three cameras are then used for steromatching to obtain 3D (three-dimensional) positions and tracks.
Camera syncing is checked by calculating reprojection errors for the three cameras. To estimate tracking error the 2D (two-dimensional) positions are stereomatched to obtain 3D positions and the 3D positions are again converted to their original 2D image plane positions to obtain a reprojection error, defined as the Figure 2. Images from the three GoPro cameras used in our experiments. The raw image as seen by the three views and a zoomed in picture of a detected honeybee with its binarized image shown below. The center of the honeybee is found using the center of mass of the bright pixels (red crosses) which is then used for steromatching and finding the 3D coordinates of the honeybee. The reprojection from the stereomatched 3D position back to the 2D camera plane is also shown to demonstrate the tracking quality (black circles). mean distance between the position of original centers found and reprojected centers across all three cameras. The frames are then adjusted such that the reprojection error is minimized and we achieve optimal camera syncing not achievable solely with the remote triggering. Furthermore, we discard any points with reprojection errors larger than 10 pixels as bad matches. In this way we are able to reduce tracking errors where for all available data, a mean imaging reprojection error of 1.6 pixels is achieved. This error is an acceptable value since each honeybee is approximately ∼10 3 pixels in size in each image and the reprojection error in finding honeybee positions can be attributed to the center of mass position identification method where the positions identified in the 2D image of the three cameras may not be exactly the same position in real-life coordinates. Additional errors may be caused by the cameras still being slightly out of sync even after optimization. The imaging and reprojection for an example frame is shown in figure 2 after optimization. Future work can improve imaging by fitting a model honeybee to the images to find the same center on each camera, which would require us to use higher resolution cameras capable of exact syncing.
After obtaining honeybee 3D positions, honeybee velocities v and accelerations a are obtained by fitting a second order polynomial to the 3D position data. The polynomial is of the form p(x) = p 1 x 2 + p 2 x + p 3 where, p 2 = v and p 1 = a. Here we only consider tracks for further analysis that are longer than 15 frames. The fit-length used for fitting polynomials is h = 7 where a track with length l is split into l − 2h sections for polynomial fitting. The rest of the analysis follows from these values.

Figures 3(a)-(c)
show the magnitude of the velocity, acceleration, and flight power, of honeybees as violin plots for when fans are off and no wind is generated and when fans are on and generating turbulent wind with the active grid. The flight power per insect mass is calculated using a formulation know from Lagrangian turbulence to define the rate of change of the kinetic energy following a tracer particle [39], where the power per mass is, P = a · v. Here we assume a negligible variance in the mass between individual honeybees compared to the large variance in acceleration and velocity. We generated wind with different turbulence levels represented by the Taylor microscale Reynolds number R λ ∼ √ Re, where Re is the Reynolds number defined as the ratio between fluid inertial forces and viscous forces. No clear trend relating the mean values of v, a, and P to the turbulence intensities produced in our experiments can be discerned in the smaller sub-plots, or from the shape of the probability distributions shown in the violin plots.
This is in agreement with previous numerical studies of bee flight where aerodynamic quantities averaged over wing-beat cycles do not change significantly in strong turbulence [28], and field experiments showing that mean bee foraging activity is not influenced by windy conditions [24].  fans are active and wind is generated and when fans are off and no generated wind is present. To better determine whether the strong directional generated wind has any effect on honeybee flight, we have combined all data with fans active into one data-set and labeled it as all wind data. We can then compare this dataset with data that was obtained with no fans on, labeled no wind data. Figures 4(a) and (b) show the probability distribution function (PDF) of honeybee velocities and accelerations. The honeybee velocity follows a normal distribution as expected. The closely normal distribution of the acceleration PDF shows that honeybees are able to maintain control while flying in turbulent conditions. Honeybees are active flyers since they are continuously flapping their wings during flight as compared to a passive flyer such as the soaring golden eagle that uses turbulence to its advantage and has acceleration PDFs more closely resembling tracers in turbulence [5], with much longer tails as compared to the honeybee. However, the accelerations while mostly close to a normal distribution, show signs of influence from the turbulent flow in the long tail of the acceleration PDF similar to particles in turbulence that reveal the intermittent nature of turbulence [2,3].
To further understand the flight dynamics, we investigate the cosine of the angle between a and v. A small angle between a and v indicates that the honeybee is acting to speed up the motion in the current direction of flight, whereas a larger angle indicates a change in flight direction. Figure 5 shows a · v/| a|| v| for all the data with no wind in (a) and all wind data in (b). From these PDFs we see that the honeybees spend the majority of their flight time in an accelerating mode, meaning that the honeybees require more time to accelerate to a certain speed, while they can rapidly decelerate and lower their velocity. This is further demonstrated by conditioning a · v/| a|| v| on the honeybee accelerations where we see that for extreme values of the acceleration, a > 2M(a) for no wind and a > 3M(a) for with wind, this value is more likely negative and that the honeybee is mostly decelerating. Here M represents the median and the a > 3M(a) data is not shown for the no wind data due to unconverged statistics in no wind conditions. Additionally, due to the unconverged statistics the bin size was adjusted accordingly for better representation of the data in no wind conditions for the a > M(a) and a > 2M(a) conditions. Even so, with the condition of a > 2M(a), the increase in deceleration likelihood is observed. This lower acceleration threshold for observing deceleration's in no wind conditions can be due to the lower power needed for the honeybee to decelerate with no wind present. This conditioning shows that deceleration events are sudden and result in large values of the acceleration.
The larger deceleration experienced by honeybees is also visible in the longer negative tails observed in the power PDF of figure 6. Here the honeybee flight power is positively skewed further showing that honeybees are accelerating most of the time, with sharp deceleration events, visible from the longer negative tail of the power PDF. The third standardized moment, or skewness, of the flight power for no wind and windy conditions are 1.74 and 3.28 respectively.
where v i is the standard deviation of v i . The acceleration a n is similarly defined.
An example of a honeybee track is shown in figure 5(c) where the slow accelerating and sudden decelerating regions are shown. Here the color shows the kinematic power of the honeybee flight. The accelerating section of the trajectory is a region where flight power slowly builds up and the honeybee increases its velocity. The decelerating section involves an abrupt turn where the power suddenly drops and the honeybee reduces its velocity as the turn is initiated followed by a power gain as the honeybee builds-up velocity once more. This phenomenon is a well-known feature for neutrally buoyant tracer particles in a turbulent flow, which experience so-called flight-crash events of strong intermittent power loss [39]. For the case of a honeybee crash events occur when abrupt turns are made.
We further examine honeybee flight by splitting tracks into flight away and towards the hive. This allows us to better analyze the influence of head-on flight against the wind (head wind) and flight in the same direction as the wind (tail wind). We determine flight direction using the direction of the honeybee velocity vector along the y coordinate axis. y is both the direction of travel away and towards the hive and the direction of the wind mean flow. Figure 7 shows a · v/| a|| v| for the split tracks for no wind data flying away from the hive (a) and towards the hive (b), and for all wind data flying away from the hive or flying against the wind (c) and flying towards the hive or flying in the same direction as the wind (d).
Interestingly, the honeybee is increasing its overall cruising speed ( a · v > 0) more frequently than decreasing it ( a · v < 0) as it leaves and flies away from the hive in both no wind and windy conditions. And for times when flying towards the hive, the accelerating and decelerating events are more evenly distributed throughout its flight, for both no wind and windy conditions, signifying its anticipation of landing back at the hive and the need to reduce velocity. Of particular interest is the comparison between figures 7(a) and (c), here the honeybee is mostly accelerating, however, in no wind conditions we clearly see two maxima when a · v/| a|| v| > 0. This indicates that when wind is not present the honeybee flies straight and makes abrupt turns further demonstrating its flight-crash like behavior [39]. Straight flight is shown by the maxima at a · v/| a|| v| > √ 2/2, and the abrupt turns are shown by the maxima at 0 < a · v/| a|| v| < √ 2/2. Here the dashed lines at ± √ 2/2 show regions where the translational acceleration a t and centripetal acceleration a c are equal, where a = a t + a c . In the case of figure 7(a), when a · v/| a|| v| > √ 2/2, the honeybee acceleration is predominantly in the form of translational acceleration a t > a c , which signifies straight flight. While the second maxima at 0 < a · v/| a|| v| < √ 2/2 shows that the honeybee is making turns and most acceleration is in the form of centripetal acceleration a c > a t . However, in windy conditions the honeybee no longer shows two distinct maxima while flying away from the hive as seen by figure 7(c) and the flight is spread across a range of a · v/| a|| v| > 0. This shows that the honeybee is flying with a range of different centripetal accelerations, where the flight trajectory is perturbed and the bee is twisting and turning with increased crash events as it flies against the wind, while translational acceleration and more straight flight is suppressed.
To quantify the increase in maneuverability of the honeybee in turbulent conditions we use the trajectory curvature k used to quantify the curvature of Lagrangian trajectories in turbulence [40]. The  curvature has also previously been used to quantify the maneuverability of midges [41]. The curvature of a trajectory is defined as k = a n /| v| 2 , where a n = | a × v|/| v| is the magnitude of the normal acceleration to the flight direction. The value for the curvature is large whenever abrupt turns are performed and such events are distinct from the background [41]. Similar to methods used by Puckett et al [41] to define significant changes in motion, we define significant changes to be instances when the curvature exceeds the mean curvature for a particular experiment. Here, the mean curvature in calm conditions is 0.025 mm −1 and for all wind data the mean curvature is 0.022 mm −1 . To determine the maneuverability change between calm and turbulent conditions, we calculate the total length of high curvature sections of trajectories in a given dataset and normalize by the total length of that dataset. In this way, we find that in calm conditions the normalized high curvature sections of motion is 0.1483 and in turbulent conditions this value is 0.1914, clearly showing the increase in honeybee maneuverability in turbulent conditions. Time-dependent statistics have also been considered where the auto-correlation function of honeybee velocity components further show the increase in turning maneuvers in windy conditions (see appendix for more details).
The PDF of a · v/| a|| v| can be conditioned on the radius of flight and is shown in figure 8 for no wind conditions and windy conditions while flying away and towards the hive. We calculate the radius of honeybee flight by obtaining the centripetal acceleration and using the relation a c = v 2 /r. Honeybee flight is not significantly affected by radii smaller than the median radius where statistics are almost unchanged for r < M(r) and r < 1/2 M(r). Thus, honeybees are most likely to fly in radii equal to or larger than that of the median radius. Here for no wind present, the median radius while flying away from the hive is 15.7 cm and while flying towards the hive is 13.8 cm. And for windy conditions the median radius while flying away from the hive is 21.9 cm and for flying towards the hive is 12.7 cm. As can be seen from figures 8(a) and (c), the maxima at smaller values of the cosine at 0 < a · v/| a|| v| < √ 2/2 in no wind conditions, as compared to windy conditions, show that the honeybee is able to make smaller radius turns when wind is not present. While, when the honeybee moves in windy conditions and flies against the wind, due to less flight control, it flies in larger radii. However, the median radius does not significantly change as the honeybee flies towards the hive in both no wind and windy conditions. During flight away from the hive in figures 8(a) and (c), we see a similar flight pattern for small radii as before. In no wind conditions during acceleration, the honeybee acceleration is either translational or centripetal. While, in windy conditions during acceleration, the translational acceleration is suppressed and accelerations are mostly centripetal due to increased turning and maneuvering. This is also observed in flight towards the hive in figures 8(b) and (d). Here, in no wind conditions, honeybee flight is evenly distributed between translational and centripetal accelerations during deceleration, while in windy conditions the honeybee has mostly centripetal acceleration during deceleration, as can be seen by the maxima at − √ 2/2 < a · v/| a|| v| < 0. Overall, in windy conditions the honeybee flight mostly consists of turning and performing evasive maneuvers for both accelerating and decelerating events during flight away and towards the hive, while in no wind conditions part of the flight also consists of more straight flight with larger contributions to translational acceleration.

Conclusions
Statistically, the honeybee velocities, accelerations, and overall flight performance are not dependent on the turbulence intensity based on the experiments conducted in this study and the turbulence intensities generated. However, further work is needed to study honeybee flight performance at larger mean wind velocities to determine how flight is affected. To observe the influence of turbulence on honeybee flight we look at the angle between the velocity and acceleration vector, the cosine of this angle reveals interesting behavior in honeybee flight. Form the cosine of this angle we see that honeybees in general spend the majority of their flight time accelerating, while deceleration events are sudden with extreme values of acceleration. This means that accelerations to top speeds take a longer time than decelerating to lower speeds, analogous to flight-crash events in turbulence [39].
The cosine of the angle between a and v further shows that in calm conditions with no wind, a honeybee is either flying straight or making turns during flight. However, as wind is introduced, straight flight is suppressed and flight is perturbed in a way that it mostly consists of turning and performing evasive maneuvers.
Future work would also need to examine statistics in the reference frame of the honeybee to investigate the effects of body motion and rotation of the honeybee in turbulence. This can be achieved with more advanced cameras and imaging at a higher frame rate and pixel density, which was not reliably possible with the accessible GoPro system used in these experiments. The results presented here on how honeybees navigate turbulence along with future studies expanding on this work can help advance miniature robotics for more stable and effective flight resistant to damage and flight path deviations due to turbulent wind. On the other hand the x component of the velocity has the shortest correlation time for both no wind and windy conditions in accordance with the unstable characteristic of roll for insect flight [6,26]. And in our experiment as seen from figure 1, the x direction is perpendicular to the main flight path which corresponds to the rolling motion of the honeybees. The y component of the velocity, v y crosses zero at an earlier time in windy conditions as compared to no wind conditions. This further demonstrates that in windy conditions honeybee flight is perturbed by increased turning maneuvers.