We present a systematic and quantitative model of huddling penguins. In this mathematical model, each individual penguin in the huddle seeks only to reduce its own heat loss. Consequently, penguins on the boundary of the huddle that are most exposed to the wind move downwind to more sheltered locations along the boundary. In contrast, penguins in the interior of the huddle neither have the space to move nor experience a significant heat loss, and they therefore remain stationary. Through these individual movements, the entire huddle experiences a robust cumulative effect that we identify, describe, and quantify. This mathematical model requires a calculation of the wind flowing around the huddle and of the resulting temperature distribution. Both of these must be recomputed each time an individual penguin moves since the huddle shape changes. Using our simulation results, we find that the key parameters affecting the huddle dynamics are the number of penguins in the huddle, the wind strength, and the amount of uncertainty in the movement of the penguins. Moreover, we find that the lone assumption of individual penguins minimizing their own heat loss results in all penguins having approximately equal access to the warmth of the huddle.

To avoid prohibitively large computations and to allow for concise analysis, our model does not account for all possible details and scenarios. Rather, our goal is to provide a simple model, based on reasonable and well defined assumptions on the geometry of the huddle and the fluid mechanics of the wind. Our model recovers important features of actual huddles such as their overall shape, downwind motion, and an equal distribution of access to the benefits of the huddle among penguins. We describe the overall framework of our model in Method, including our assumptions. In Results and Discussion, we show simulation results, and identify, describe, and quantify the dynamics of the huddle. We discuss how these results compare to field observations and outline how one may extend this mathematical model to include a number of particular effects in Conclusion.

We introduce here a systematic and quantitative mathematical model for penguin huddles. This mathematical model is aligned with the qualitative observations by Le Maho stated above. Moreover, it is consistent with the idea that penguins huddle tightly to reduce their cold-exposed body surfaces, and increase the ambient temperature. The key assumption of our mathematical model is that each individual penguin seeks to reduce its own heat loss. Thus, a penguin on the boundary of the huddle exposed to the wind will move downwind along the huddle boundary. In contrast, the penguins in the interior of the huddle neither have space to move nor experience significant heat loss, so they remain stationary. While penguins inside the huddle have been observed to make multiple small displacements [13] , we consider here that these motions are small compared to the motion observed on the edge of the huddle. The accumulation of individual penguin movements along the boundary leads to coherent and robust huddle dynamics that we identify, describe, and quantify. In particular, we find from our simulation results that the number of penguins in the huddle, the wind strength quantified by the Péclet number, and the amount of uncertainty in the penguin movement are the key factors governing the dynamics of a huddle. Our simulation results show that the features of this model are sufficient to result in each penguin having equal access to the benefits of the huddle.

The huddles are not motionless; movement is extremely slow, but continuous. The huddle is urged along by the wind, the rear-flank birds (those most exposed to the wind) advancing slowly along the sides of the huddle in order to be protected from the wind. Thus, birds that at first are in the center of the huddle become members of the rear flank and move, in their turn, up the sidelines.

An important feature of huddles is that each penguin has approximately equal opportunity to the warmth of the huddle. How each penguin obtains this equal access is thought to be the result of a complex phenomenon in which penguins reorganize themselves within the huddle [1] – [3] . Gilbert et al [1] attribute heterogeneity of the huddle shape to ensuring this equal access, but without providing details as to how equality is achieved. On the other hand, Zitterbart et al [13] use ideas from condensed matter physics to explain how penguins within a huddle reorganize themselves. An important set of observations regarding penguin movement within huddles reported by Le Maho [5] is

Measurements of the body temperature of penguins in various environments and their relation to weather conditions have led to significant insight into huddle formation. [1] , [3] , [7] Huddles occur more frequently at lower ambient temperatures and in higher wind speed, but the intensity of the huddle (i.e. the number density of the huddle) depends on lower ambient temperatures only [3] . Few theoretical models of huddling have been presented. Some modeling effort has assumed that penguins move from the windward to the leeward side of the huddle, but without providing much justification [11] . Canals and Bozinovic [12] modeled huddle formation in mice (Mus musculus) as a self-organizing event. In particular, they modeled huddling as a second-order phase transition triggered by cold temperatures.

Emperor penguins huddle to conserve energy, which is particularly important since they must fast for periods of 105 to 115 days [6] . Emperor penguins benefit from huddles because of the reduction of body surface area exposed to the cold and owing to the warm temperature inside the huddle [9] . The ambient temperature in the huddle is at least 20°C and may reach as high as 37.5°C [1] . Despite the fact that huddles achieve these high ambient temperatures, emperor penguins benefit most from the huddle through the reduction of cold-exposed body surfaces [10] .

Emperor penguins (Aptenodytes forsteri) are known to huddle to survive long periods of fast in the severe conditions of Antartcic winters. They are able to form huddles because they are not tied to a fixed nest. Huddles are discontinuous events that last for relatively short durations (on the order of a few hours) [1] corresponding to storm events [2] , [3] . The number density in a huddle at a colony may be as high as m 2 [4] , [5] . Researchers have observed directly penguins huddling at their colonies [1] , [2] , [5] – [7] , and there is evidence indicating that emperor penguins also huddle during foraging trips [8] .

Methods

In our mathematical model, individual heat preservation is the goal of each and every penguin in the huddle, irrespective of the heat preservation of the huddle as a whole. In other words, each penguin will move or stay stationary to minimize its own heat loss. This mathematical model therefore requires the determination of the local rate of heat loss for each penguin in the huddle. This local rate of heat loss depends on the temperature distribution outside the huddle which, in turn, depends on the wind flow around the huddle.

Specifically, we focus on the dynamics of a single huddle, where all the penguins present are part of the huddle. This huddle is assumed to be situated on a flat plane, so there are no obstacles impeding penguin movement. Penguins in this huddle have uniform size and shape. Our model does not account for all heat exchanges between penguins and their environments. In particular, penguins are known to lose a large quantity of heat through their feet and eyes [14]. In addition, the flow of wind over the top of the huddle contributes to cooling penguins. However, these heat loses affect penguins equally, irrespective of their position within a huddle, and therefore do not have a strong influence on the huddle dynamics. In contrast, the wind flow around the sides of the huddle affects penguins differently depending on whether they are on the edge or near the center of the huddle. For this reason, we model only the wind flowing on a two dimensional plane around the area occupied by the huddle. Our procedure to simulate huddles is as follows:

Generate a huddle and determine the huddle boundary. Compute the wind flow around the huddle. Compute the temperature profile around the huddle. Compute the local rate of heat loss for each penguin. Add random variations to the rate of heat loss (optional). Identify the penguin with the highest rate of heat loss (the “mover”) and move it to a location on the boundary where heat loss is minimal. Determine the new huddle boundary. Repeat over the desired number of iterations by going back to Step 2.

In what follows, we describe in detail each step of this procedure. Included in this discussion are any simplifying assumptions, and their justification.

1. Generate a Huddle and Determine the Huddle Boundary From observing videos of huddling penguins [13], [15], we note that huddling penguins are packed very tightly, presumably to minimize cold-body surfaces and maximize the ambient temperature in the huddle. It is well known mathematically that the densest packing of disks on a flat plane corresponds to having the center of each disk placed on a hexagonal or honeycomb lattice, so that each disk is inscribed in a regular hexagon and touches six neighboring disks. Moreover, direct observations reveal that huddling penguins generally position themselves on a hexagonal grid [13]. For these reasons, we assume in our model that penguins in the huddle stand centered on a hexagonal grid. Moreover, all points on the hexagonal lattice corresponding to the huddle’s interior are assumed to be occupied, so that there are no empty spaces within the huddle. Thus, the huddle boundary is determined uniquely by connecting the locations of penguins with fewer than six neighbors. Furthermore, we assume that all penguins in the huddle have at least two neighbors, in which case, the area generated by connecting the lattice points on the huddle boundary is a polygon. We therefore initiate our simulations by generating a huddle satisfying these conditions, starting with five penguins and adding penguins at locations chosen randomly among the eligible positions (adjacent to the huddle, with two neighbors, and leaving no empty space). Once the initial huddle is formed, we consider that the number of penguin within it remains constant.

2. Compute the Wind Flow Around the Huddle To determine the wind flow around the huddle, we need only to consider a two dimensional flow around a polygon. Moreover, we assume that this flow is inviscid and irrotational. These assumptions imply that we do not resolve turbulent wind flows, which would obfuscate the computation of the temperature profile needed to compute the local rate of heat loss. Rather, we find a smooth, regular wind flow around the huddle. Nonetheless, the key relationships in this mathematical model between the wind, the temperature, and individual heat loss are not compromised by these assumptions. Because the wind flow is significantly faster that the movement of a penguin, we assume that this flow is steady. In other words, the wind flow does not depend on the time elapsed since a penguin has relocated, and is only dependent on the huddle shape. Consequently, we are able to use the mathematics of complex variables and the physical theory of potential flow to describe the flow around the huddle. Let denote the wind velocity potential, and the wind velocity itself, , be the gradient of : . We combine the corresponding streamlines, denoted by , with to form the complex potential, defined as with i denoting the imaginary constant. Because the flow is irrotational, F is an analytic function on the flat plane outside the huddle. Since the huddle boundary is a polygon, we can make use of the Schwarz-Christoffel transformation, denoted by . The Schwarz-Christoffel transformation is a conformal mapping from the outside of a disk, in what we will call the canonical domain denoted by , to the outside of a polygon, in what corresponds to the physical domain, denoted by [16]. In the context of our model, the polygon in the physical domain corresponds to the boundary of the penguin huddle. This transformation is useful here because in the canonical domain, we are able to compute F easily via the Joukowsky transform , where U is the wind speed away from the disk. This transform maps the portion of the upper half-plane of the canonical domain outside the unit circle to the upper-half plane corresponding to with denoting the imaginary part. In this transformed domain, the complex potential is known exactly. Mathematically, the wind flow around the polygon is therefore given by with denoting the real part of a complex valued function and is the inverse of the Schwarz-Christoffel transformation. The key step in this computation is the Schwarz-Christoffel transformation. Because huddle shapes may be arbitrary polygons, we are not generally able to compute analytically. Instead, we use the freely available software package developed by Driscoll [17] which computes this transformation numerically. We then obtain a description of the wind flow at any point outside the huddle. Note that the wind speed, , only affects the magnitude of the flow, but does not modify the location of the streamlines.

3. Compute the Temperature Profile Around the Huddle The temperature far outside the huddle, denoted by , is significantly less than the temperature inside the huddle interior. The large transition in temperature from outside the huddle to inside the huddle is most significant near the huddle boundary where penguins are most exposed to the wind. To model this situation in our simulations, we assume that the temperature on the edge of the huddle, , and are constants with . This assumption allows us to focus our attention solely on the temperature profile in the transition region near the huddle boundary. Using the wind velocity around the huddle, , we solve for the temperature profile obtained from the advection-diffusion equation. As was the case for the wind flow, we assume that the temperature profile is steady, meaning that the temperature profile reaches equilibrium faster than any penguin motion. Consequently, the spatial distribution of the temperature around the huddle satisfies the steady advection-diffusion equation where is the diffusivity of the temperature. This equation is invariant under conformal mapping provided may be written as the gradient of a potential function. We can therefore use the Schwarz-Christoffel mapping once again and solve for the temperature around a disk in the canonical domain, before mapping it back around the polygonal huddle boundary. To simplify the analysis, we non-dimensionalize the equations by rescaling the velocity by the wind velocity away from the huddle , as . We rescale the temperature by the temperature difference between the penguins and the air away form the huddle , so that . We also rescale distances by the size of the huddle , which corresponds to the radius of a huddle with the same number of penguins disposed to fill a disk, as , . Rewriting the equation above in terms of these new variables, we obtain the non-dimensional equation with denoting the Péclet number. This non-dimensional number is proportional to the wind speed and therefore captures the effects of varying the strength of the wind. The boundary conditions we use to solve this equation in the canonical domain are on the unit circle, where the penguins temperature remains effectively constant, and on a circle of radius , away from the huddle. We compute the temperature distribution exterior to the unit circle in the canonical domain by discretizing the annulus between the unit circle and the circle of radius using finite differences. In particular, we discretize the annulus into regions of equal radius and angular increments. We then replace the differential operators in the equation with centered, finite difference approximations. The result leads to a linear system that we solve numerically. Using the Schwarz-Christoffel mapping, we may then recover the temperature distribution around the huddle. As an example, we show in Figure 1 a sample of the temperature profile obtained around a huddle made up of 100 penguins, choosing . PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 1. Sample of a computed temperature. Temperature distribution around a huddle of 100 penguins, for . Here red and blue correspond to warmer and cooler temperatures, respectively. Individual penguins are shown in black, as is the boundary of the huddle, while the polygonal interior of the huddle is shown in white. We show an initial configuration where penguins have yet to be relocated. https://doi.org/10.1371/journal.pone.0050277.g001 When the flow around the huddle is turbulent, as is likely to be the case for all but the mildest winds [18], the relevant diffusivity is a turbulent diffusivity. For that case, the relevant Péclet number is therefore a turbulent Péclet number. While this quantity varies depending on the flow, one may estimate it as being analogous to a turbulent Reynolds number, , which is based on the diffusivity of momentum rather than temperature. Because the diffusion in turbulent flow is due to correlated fluid motions that mix temperature and momentum in similar manners, the turbulent Reynolds and Péclet numbers have similar values. The turbulent Reynolds number may be estimated as the value where the drag coefficient becomes effectively constant, which typically occurs for [19]. In the present study, we therefore limit our attention to flows with .

4. Compute the Local Rate of Heat Loss for Each Penguin Penguins at the huddle boundary experience the most significant heat loss compared to those within the huddle interior since the temperature profile changes most abruptly outside the huddle boundary. We seek to find the penguin on the huddle’s edge with the largest rate of heat loss. Therefore, we compute the local rate of heat loss only for penguins on the huddle boundary. The local rate of heat loss at a boundary is proportional to the derivative of the temperature in the direction normal to the boundary [20]. From our computed temperature profile, we may approximate the normal derivative along the unit circle (parameterized by angle ) in the canonical domain using where is the distance between consecutive points in the radial direction. To determine the heat loss associated to each penguin on the boundary, we first define points that are located halfway between consecutive penguins on the boundary, respectively labeled as penguins and . We then find the preimage of these points along the unit circle as . The heat loss associated to penguin is then given by We approximate using the numerically obtained temperature profile obtained as described above to approximate the normal derivative, which we then integrate numerically. We note that is always positive, indicating that heat is always lost. Penguins inside the huddle are assumed to lose only a negligible amount of heat compared to those on the periphery of the huddle, and we therefore set for penguins that are not on the edge of the huddle.

5. Add Random Variations (Optional) The heat loss computation described above is idealized because everything is assumed to be known with absolute certainty. To account for variations that can occur in real huddles, we add uncertainty to this model through a random perturbation to the heat loss associated to each penguin on the boundary. We let denote the average heat loss experienced by penguins disposed in a disk subject to the same conditions as the model huddle, with the number of penguins in the huddle. We introduce uncertainty in our mathematical model by assigning to each penguin an effective heat loss that includes a random component where is a random number drawn from a uniform distribution ranging from −0.5 to 0.5, and is a parameter quantifying the magnitude of the random component relative to that of the heat lost to the ambient air. While remains constant throughout the iterative process, is re-drawn each time a penguin moves and the heat loss is recomputed. We then use the effective heat loss, , to determine the dynamics of the huddle. The parameter determines the importance of random variations in the system. The special case yields a completely deterministic, and idealized, system. No individual variations are captured in this case, but the results are perfectly replicable, which simplifies the analysis. As increases, the effects of random variations increase, and the importance of wind-related heat loss progressively diminishes. For large values of , say for , the wind direction would become negligible and model penguins would perceive a heat loss that is essentially random. It is likely that the relevant magnitude of depends on the harshness of the conditions to which the huddle is subjected. In low winds or relatively warm temperature, the heat lost by each penguin is not very large, and penguin behavior is therefore likely to contain more random variations. On the other hand, in very harsh conditions, the imperative to minimize heat loss will be much stronger, and departure from the optimal behavior, in the form of random variations, are less likely.