Access the full text.
Sign up today, get DeepDyve free for 14 days.
Richard Gerum, Ben Fabry, Claus Metzner, Michaël Beaulieu, A. Ancel, D. Zitterbart (2013)
The origin of traveling waves in an emperor penguin huddleNew Journal of Physics, 15
D. Crowdy, Hyun-Joon Kang (2001)
Squeeze flow of multiply-connected fluid domains in a Hele-Shaw cellJournal of Nonlinear Science, 11
W.W. Mullins (1964)
10.1063/1.1713333J. Appl. Phys., 35
Pablo Brubeck, Y. Nakatsukasa, L. Trefethen (2019)
Vandermonde with ArnoldiArXiv, abs/1911.09988
J. Ockendon (2004)
The classical Stefan Problem: Basic Concepts, Modelling and Analysis. By S. C. GUPTA. Elsevier, 2003. 385 pp. ISBN 0444 510869. £66.50Journal of Fluid Mechanics, 520
F. Dutka, V. Starchenko, F. Osselin, S. Magni, P. Szymczak, A. Ladd (2020)
Time-dependent shapes of a dissolving mineral grain: Comparisons of simulations with microfluidic experimentsChemical Geology
(1986)
SD851674 Howison
J. Hilton, Claire Miller, J. Sharples, A. Sullivan (2016)
Curvature effects in the dynamic propagation of wildfiresInternational Journal of Wildland Fire, 25
C. Gilbert (2006)
10.1016/j.physbeh.2006.04.024Physiol. Behav., 88
S. Richter, Richard Gerum, A. Winterl, A. Houstin, M. Seifert, J. Peschel, B. Fabry, C. Bohec, D. Zitterbart (2018)
Phase transitions in huddling emperor penguinsJournal of Physics D: Applied Physics, 51
C. Herreid (1963)
Temperature Regulation of Mexican Free-Tailed Bats in Cave HabitatsJournal of Mammalogy, 44
Tamzidul Mina, B. Min (2018)
Penguin Huddling inspired Distributed Boundary Movement for Group Survival in Multi-robot Systems using Gaussian Processes2018 IEEE International Conference on Robotics and Biomimetics (ROBIO)
C.L. Williams, J.C. Hagelin, G.L. Kooyman (2015)
Hidden keys to survival: the type, density, pattern and functional role of emperor penguin body feathersProc. R. Soc. Lond. B, 282
S. Costa, L. Trefethen (2021)
AAA-least squares rational approximation and solution of Laplace problemsArXiv, abs/2107.01574
A. Ancel, C. Gilbert, N. Poulin, M. Beaulieu, B. Thierry (2015)
New insights into the huddling dynamics of emperor penguinsAnimal Behaviour, 110
M. Bazant, D. Crowdy (2004)
Conformal mapping methods for interfacial dynamics
D. Zitterbart, B. Wienecke, J. Butler, B. Fabry (2011)
Coordinated Movements Prevent Jamming in an Emperor Penguin HuddlePLoS ONE, 6
L. Trefethen (2019)
Numerical Conformal Mapping with Rational FunctionsComputational Methods and Function Theory, 20
Aaron Waters, F. Blanchette, A. Kim (2012)
Modeling Huddling PenguinsPLoS ONE, 7
Wen Gu, Jason Christian, C. Woodson (2018)
A novel coupled fluid-behavior model for simulating dynamic huddle formationPLoS ONE, 13
C. Dawson, J. Vincent, G. Jeronimidis, G. Rice, Paula Forshaw (1999)
Heat transfer through penguin feathersJournal of theoretical biology, 199 3
C. Ryan, Lynne Burns, H. Broders (2019)
Changes in underground roosting patterns to optimize energy conservation in hibernating batsCanadian Journal of Zoology
J. Sethian (1985)
Curvature and the evolution of frontsCommunications in Mathematical Physics, 101
Peter Baddoo (2020)
Lightning Solvers for Potential FlowsFluids
S. Bernardi, M. Scianna (2020)
An agent-based approach for modelling collective dynamics in animal groups distinguishing individual speed and orientationPhilos. Trans. R. Soc. Lond. B, 375
L. Paterson (1981)
Radial fingering in a Hele Shaw cellJournal of Fluid Mechanics, 113
C. Gilbert, G. Robertson, Y. Maho, Y. Naito, A. Ancel (2006)
Huddling behavior in emperor penguins: Dynamics of huddlingPhysiology & Behavior, 88
S. Howison (1986)
Fingering in Hele-Shaw cellsJournal of Fluid Mechanics, 167
Linda Cummings, Y. Hohlov, S. Howison, K. Kornev (1999)
Two-dimensional solidification and melting in potential flowsJournal of Fluid Mechanics, 378
The two-phase melting/freezing problem: exterior and interior effects. Porous media ﬂow about freeze pipes where interior and exterior temperature gradients are in effect
O. Agam (2008)
Viscous fingering in volatile thin films.Physical review. E, Statistical, nonlinear, and soft matter physics, 79 2 Pt 1
G. Kooyman, R. Gentry, W. Bergman, H. Hammel (1976)
Heat loss in penguins during immersion and compression.Comparative biochemistry and physiology. A, Comparative physiology, 54 1
V. Tsai, J. Wettlaufer (2007)
Star patterns on lake ice.Physical review. E, Statistical, nonlinear, and soft matter physics, 75 6 Pt 2
Gary Nave, Nelson Mitchell, Jordan Dick, Tyler Schuessler, Joaquin Lagarrigue, O. Peleg (2019)
Attraction, Dynamics, and Phase Transitions in Fire Ant Tower-BuildingFrontiers in Robotics and AI, 7
M. Dallaston, S. McCue (2016)
A curve shortening flow rule for closed embedded plane curves with a prescribed rate of change in enclosed areaProceedings. Mathematical, Physical, and Engineering Sciences / The Royal Society, 472
Y. Nakatsukasa, O. Sète, L. Trefethen (2016)
The AAA Algorithm for Rational ApproximationSIAM J. Sci. Comput., 40
Cassondra Williams, J. Hagelin, G. Kooyman (2015)
Hidden keys to survival: the type, density, pattern and functional role of emperor penguin body feathersProceedings of the Royal Society B: Biological Sciences, 282
P. Grodzki, P. Szymczak (2019)
Reactive-infiltration instability in radial geometry: From dissolution fingers to star patterns.Physical review. E, 100 3-1
The “full" penguin problem: exterior, interior and area conservation effects. This is the full huddle continuum problem covered in this work
P. Basu (2013)
Biomass Gasification, Pyrolysis and Torrefaction: Practical Design and Theory
T. Driscoll (1996)
Algorithm 756: a MATLAB toolbox for Schwarz-Christoffel mappingACM Trans. Math. Softw., 22
L. Trefethen (2018)
SERIES SOLUTION OF LAPLACE PROBLEMSThe ANZIAM Journal, 60
R. Kirkwood, G. Robertson (1999)
The Occurrence and Purpose of Huddling by Emperor Penguins During Foraging TripsEmu, 99
G. Taylor, P. Saffman (1959)
A NOTE ON THE MOTION OF BUBBLES IN A HELE-SHAW CELL AND POROUS MEDIUMQuarterly Journal of Mechanics and Applied Mathematics, 12
Sara Bernardi, M. Scianna (2020)
An agent-based approach for modelling collective dynamics in animal groups distinguishing individual speed and orientationPhilosophical Transactions of the Royal Society B, 375
Richard Gerum, S. Richter, B. Fabry, C. Bohec, F. Bonadonna, A. Nesterova, D. Zitterbart (2018)
Structural organisation and dynamics in king penguin coloniesJournal of Physics D: Applied Physics, 51
S. Harris, N. McDonald (2022)
Fingering instability in wildfire frontsJournal of Fluid Mechanics, 943
A. Gopal, L. Trefethen (2019)
Solving Laplace Problems with Corner Singularities via Rational FunctionsSIAM J. Numer. Anal., 57
S. Costa (2020)
Solving Laplace problems with the AAA algorithmArXiv, abs/2001.09439
W. Mullins (1964)
Stability of a Planar Interface During Solidification of a Dilute Binary Alloy
P. Saffman, G. Taylor (1958)
The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquidProceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 245
Hungtang Ko, Ting-Ying Yu, D. Hu (2022)
Fire ant rafts elongate under fluid flowsBioinspiration & Biomimetics, 17
A. Ladd, Liang Yu, P. Szymczak (2020)
Dissolution of a cylindrical disk in Hele-Shaw flow: a conformal-mapping approachJournal of Fluid Mechanics, 903
M. Burger, J. Haskovec, Marie-Therese Wolfram (2013)
Individual based and mean-field modeling of direct aggregationPhysica D. Nonlinear Phenomena, 260
M. Goldstein, R. Reid (1978)
Effect of fluid flow on freezing and thawing of saturated porous mediaProceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 364
Y. Le Maho (1977)
The emperor penguin: a strategy to live and breed in the cold: morphology, physiology, ecology, and behavior distinguish the polar emperor penguin from other penguin species, particularly from its close relative, the king penguinAm. Sci., 65
A Plummer (1967)
Introduction to Solid State PhysicsPhysics Bulletin, 18
C. Rycroft, M. Bazant (2015)
Asymmetric collapse by dissolution or melting in a uniform flowProceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472
R. McDonald, M. Mineev-Weinstein (2015)
Poisson growthAnalysis and Mathematical Physics, 5
Kunal Bhattacharya, T. Vicsek (2010)
Collective decision making in cohesive flocksNew Journal of Physics, 12
Jaehyuk Choi, D. Margetis, T. Squires, M. Bazant (2004)
Steady advection–diffusion around finite absorbers in two-dimensional potential flowsJournal of Fluid Mechanics, 536
J. Langer (1980)
Instabilities and pattern formation in crystal growthReviews of Modern Physics, 52
V. Entov, P. Etingof (1991)
BUBBLE CONTRACTION IN HELE-SHAW CELLSQuarterly Journal of Mechanics and Applied Mathematics, 44
R. Brower, D. Kessler, J. Koplik, H. Levine (1984)
Geometrical models of interface evolutionPhysical Review A, 29
Björn Gustafsson, A. Vasiliev (2006)
Conformal and Potential Analysis in Hele-Shaw Cells
Y. Lemaho (1977)
The Emperor Penguin: A Strategy to Live and Breed in the Cold, 65
Dominic McCafferty, Caroline Gilbert, A. Thierry, A. Thierry, John Currie, Y. Maho, Y. Maho, A. Ancel, A. Ancel (2013)
Emperor penguin body surfaces cool below air temperatureBiology Letters, 9
M. Dallaston, S. McCue (2013)
An accurate numerical scheme for the contraction of a bubble in a Hele-Shaw cell
A. Harker (2018)
The classical Stefan problem: basic concepts, modelling and analysis with quasi-analytical solutions and methodsContemporary Physics, 59
Y. Katz, K. Tunstrøm, C. Ioannou, C. Huepe, I. Couzin (2011)
Inferring the structure and dynamics of interactions in schooling fishProceedings of the National Academy of Sciences, 108
Penguins huddling in a cold wind are represented by a two-dimensional, continuum model. The huddle boundary evolves due to heat loss to the huddle exterior and through the re- organisation of penguins as they seek to regulate their heat production within the huddle. These two heat transfer mechanisms, along with area, or penguin number, conservation, gives a free boundary problem whose dynamics depend on both the dynamics interior and exterior to the huddle. Assuming the huddle shape evolves slowly compared to the advective timescale of the exterior wind, the interior temperature is governed by a Poisson equation and the exterior temperature by the steady advection-diffusion equation. The exterior, ad- vective wind velocity is the gradient of a harmonic, scalar ﬁeld. The conformal invariance of the exterior governing equations is used to convert the system to a Polubarinova-Galin type equation, with forcing depending on both the interior and exterior temperature gradi- ents at the huddle boundary. The interior Poisson equation is not conformally invariant, so the interior temperature gradient is found numerically using a combined adaptive Antoulas- Anderson and least squares algorithm. The results show that, irrespective of the starting shape, penguin huddles evolve into an egg-like steady shape. This shape is dependent on the wind strength, parameterised by the Péclet number Pe, and a parameter β which effectively measures the strength of the interior self-generation of heat by the penguins. The numerical method developed is applicable to a further ﬁve free boundary problems. Keywords Free boundary problems · Conformal mapping · AAA-least squares algorithm · Fluid mechanics · Collective behaviour · Huddling 1Introduction Among all the penguin species living in the Antarctic region [1], the emperor penguin (Aptenodytes forsteri) experiences the most extreme weather conditions, enduring wind −1 ◦ speeds over 30 ms at temperatures below −40 C[1–4]. Unlike other species such as the gentoo and the Adélie [5], emperor penguins have no ﬁxed nest [6, 7], allowing the birds to huddle together for warmth in severe wind conditions. This strategy is very effective; ambient temperatures within the huddle often exceed 20 C, with a maximum temperature of 37.5 C having been recorded [8]. Such huddling is crucial for the penguins’ survival, S.J. Harris sam.harris.16@ucl.ac.uk Department of Mathematics, University College London, Gower Street, London, WC1E 6BT, UK 7 Page 2of17 S.J.Harris, N.R. McDonald especially during foraging trips [9] and throughout their breeding season, which coincides with the Antarctic winter [3]. Over the few hours that huddling events typically last [6, 8], the penguins are contin- uously (albeit slowly) on the move. The birds constantly reorganise their position in the huddle, where [1] observes that “birds that at ﬁrst are in the center of the huddle become members of the rear ﬂank [the side most exposed to the wind] and move, in their turn, up the sideline”. This process beneﬁts all members of the huddle and allows the heat to be equally shared: those at the edge do not remain cold and those at the centre do not overheat. These individual penguin movements affect the propagation of the huddle as a whole - the huddle is seen to propagate downwind [10] and the shape traced out by the edge of the huddle, or huddle boundary, may change over time [4, 6]. It is the propagation and evolution of the huddle boundary shape that is the focus of this work. The modelling of penguin huddle dynamics has received increased attention over the last decade, but still remains sparse in the literature, despite numerous observations and ﬁndings from the ﬁeld, see [1–3, 10]. The works [4, 7, 11] draw the analogy between penguin hud- dle/colony formation and condensed matter physics, with Lennard-Jones-like forces used to represent attractive and repulsive interactions between penguins that allow complex, lattice structures (such as huddles) to form. This idea is expanded upon in [12], where Gaussian processes are used to test if robots could adopt this penguin inspired huddling strategy. Sim- ilarly, [13] liken the rearrangement of the penguin huddle to a phase transition from solid (dense huddle) to liquid (loose huddle) to gas (no huddle) depending on environmental fac- tors such as temperature and wind speed. A rigorous ﬂuid dynamics approach is taken in [14], where a ﬁnite difference method is used to solve the full Navier-Stokes equations of the surrounding ﬂuid (the exterior wind), giving a two-way coupling between the environ- ment and the behaviour of the huddle. Using a simpler ﬂuids based model, [6] assume the exterior wind to be a two-dimensional, irrotational ﬂow of an incompressible ﬂuid and from this model, the temperature proﬁle exterior to the huddle and individual penguin movements can be found. In this work, the approach of [6] is expanded upon with two key differences. First, a simple model for the effect of interior penguin reorganisation as a result of self-generation of heat throughout the entire huddle, as observed in [1], is included. In [6], while penguin reorganisation around the huddle boundary was acknowledged, it was assumed that interior penguins were so tightly packed that they could not change their position in the huddle. That assumption is here dropped, yet it is still assumed that the penguins are packed tightly enough such that wind does not permeate through the huddle and so the wind ﬂow remains entirely in the exterior region. Throughout this work, the phrases “exterior only penguin problem” and “full penguin problem” are used to refer to the penguin huddle problem with- out and with interior penguin reorganisation, respectively. Second, the penguin huddle itself is treated as a continuum. Existing huddle models treat penguins as discrete elements (or particles) such that individual penguin movements can easily be tracked. For sufﬁciently large huddles on the order of a thousand penguins [1], these discrete models quickly become computationally expensive. Continuum models, as used here, are often faster, cheaper, and preferred when studying the change in shape and position of the overall huddle. Using a continuum to model the interior and exterior of the full penguin problem draws analogues to other two-dimensional free boundary problems. Foremost among these is the classical Hele-Shaw free-boundary problem and its extensions, including Saffman-Taylor ﬁngering (a consequence of injecting a ﬂuid into another, more viscous ﬂuid) [15–17], the expansion and contraction of bubbles in a Hele-Shaw cell [18–21] and the dissolution of soluble objects in two-dimensional potential ﬂow [22]. Using conformal mapping to re- formulate the free boundary problem as a Polubarinova-Galin (PG) equation has proved Penguin Huddling: A Continuum Model Page 3 of 17 7 effective in ﬁnding exact and numerical solutions for this class of problems. Solidiﬁcation and melting problems have also been extensively studied using analogous techniques, see [23–28]and [29]; the latter of these examines the freezing of porous media ﬂow by freeze pipes, with the velocity of the free boundary dependent on both interior and exterior temper- ature gradients, similar to the full penguin problem. Free boundary problems arise in other natural processes, including, but not limited to, dendritic crystal growth [30, 31]and ﬂame front propagation [32–34]. Continua may also be used to model other collective behaviours [4, 35, 36] such as ﬁsh schooling [37], bird ﬂocking [38], bat swarming [39, 40]and ant colony formation [41, 42]. There are two heat transfer mechanisms in the full penguin problem, which formulate a suitable boundary condition for the propagation of the huddle boundary. The ﬁrst is the cooling effect of the wind on the huddle. As in [6], the incompressible wind ﬂow is assumed irrotational and ﬂowing sufﬁciently faster than the propagation of the huddle, thus, the ﬂow can be treated as steady. As such, the wind ﬂow and temperature transport in the region exterior to the penguin huddle can be modelled by the Laplace and the steady advection- diffusion equations, respectively. This coupled pair of partial differential equations (PDEs) is conformally invariant [43] - a property which motivates the use of a conformal mapping approach in the numerical method [21–24, 34]. The exterior region and the huddle bound- ary in the physical z−plane can be conformally mapped to the exterior and boundary of the unit disk in some “mathematical” (or “canonical”) ζ −plane, and vice versa. Therefore, the physical coordinates can be expressed in terms of the conformal map z = f(ζ, t) and the condition on the huddle boundary rewritten as a PG type equation [44]. The task then becomes that of ﬁnding the conformal map. The second heat transfer mechanism is that of the uniform, self-generation of heat by the penguins in the interior. Its steady diffusion amongst the densely packed penguins is modelled by a Poisson equation for the interior temperature with a constant forcing term. By the Riemann mapping theorem (see e.g. [43]), the interior of the penguin huddle will not also be mapped to the interior of the ζ −disk by the conformal map z = f(ζ, t) used to map the exteriors. Further, the Poisson equation is not conformally invariant [45]. Hence, to solve for temperature in the interior using a conformal mapping approach, a different conformal map f would need to be found numerically and the transformed Poisson equation in the ζ − plane solved. Alternatively, the interior temperature can be found directly in the z−plane without the use of conformal mapping. Here, the methods developed by Trefethen and colleagues [46–50] are used to ﬁnd a rational approximation of a harmonic function in some simply connected domain as the real part of some analytic function, subject to given boundary conditions. This function consists of a power series in z and a series of rational functions involving poles in the exterior clustered near singularity points [47]. Both the adaptive Antoulas-Anderson algorithm (AAA, pronounced triple-A) [49] and a least- squares (LS) algorithm, such as the lightning Laplace solver [48] can be used to numerically ﬁnd this analytic function [50] and hence the required harmonic function. In this work, by ﬁrst ﬁnding a particular solution, the problem of solving the interior Poisson equation then becomes one of solving the Laplace equation with a modiﬁed boundary condition. In turn, this is readily solved numerically by the AAA-LS method. The structure of this paper is as follows. The two-dimensional, continuum penguin hud- dle model is formulated in Sect. 2, with assumptions stated and dimensionless quantities introduced. A suitable boundary condition for the normal boundary velocity v is derived, incorporating effects of exterior wind cooling, internal heat generation by the penguins and area (penguin) conservation. The numerical method is then developed in Sect. 3,where the equation for v is rewritten as a PG type equation using a conformal mapping approach. n 7 Page 4of17 S.J.Harris, N.R. McDonald Attention here is given to the exterior heat ﬂux σ in the ζ −plane and the interior heat ﬂux which is found in the physical plane. Results are presented in Sect. 4, illustrating the evo- lution of penguin huddles of different starting shapes under a variety of wind and interior heating effect strengths. Conclusions are discussed in Sect. 5 and extensions to the work, and applications to further (free boundary) problems, noted. 2ModelSetup Consider a huddle of penguins with horizontal length scale L where the penguins are as- sumed homogeneous in size and shape. The huddle is considered as a continuum with pen- guins packed as tightly as possible while still allowing for huddle reorganisation [11]and hence uniform spacing between each penguin can be assumed constant across the entire huddle and for all time. Emperor penguins have a typical height of H = 1m [1]and their tight packing means that background winds are diverted around the huddle, allowing the dy- namics to be considered as two-dimensional, see Fig. 1. The boundary of the penguin huddle γ(t) is deﬁned as the two-dimensional, non-overlapping curve enclosing the huddle; the ﬁ- nite, simply connected interior of the curve occupied by the huddling penguins is labelled R and the penguin-free exterior as . The boundary curve γ is orientated anti-clockwise such that the unit normal n ˆ points from the interior to the exterior. The huddle R is buffeted by a unidirectional and constant strength wind of magnitude U which, without loss of generality, ﬂows in the positive x ˆ direction. The wind is the irrota- tional ﬂow of air, which is considered to be an inviscid and incompressible ﬂuid which ﬂows sufﬁciently rapidly compared to the evolution of the penguin huddle as a whole, so that the transport of temperature can be treated as steady [6]. In reality, viscous effects taking the form of thin boundary layers around the exterior of the penguin huddle are present. In this work, these are ignored in order to facilitate the accurate computation of the exterior po- tential ﬂow using conformal mapping based methods. Denote the temperature far from the huddle as T , the temperature on the huddle boundary as T and the maximum temperature ∞ f in the huddle interior as T . The equations governing the wind and exterior temperature transport are, respectively, the Laplace equation and the steady advection-diffusion equation ∇ φ =0in , (1) u ·∇T = D ∇ T in , (2) where φ is the wind velocity potential, u = ∇φ is the wind velocity, T is the temperature in and D is the diffusivity of temperature in air. The associated boundary and far ﬁeld Fig. 1 Schematic diagram of the full penguin problem showing the interior region R and exterior region . a) Plan view where T denotes the maximum temperature T inside the huddle, b) side view c Penguin Huddling: A Continuum Model Page 5 of 17 7 conditions are ∂φ T = T , =0on γ, (3) ∂n T → T ,φ → Ux as r →∞. (4) Now, consider the interior R of the huddle where there is no wind. Owing to their ability to self-generate heat, which subsequently diffuses within the huddle, penguins nearer the huddle centre are warmer than those at the boundary [1, 11], i.e. T >T . This source of c f heat provides a further contribution to the heat ﬂux at the boundary γ . Consistent with the continuum assumption, it is assumed that heat is generated uniformly across the huddle. Further, as in the exterior, it is assumed the interior heat diffuses quickly compared to the evolution of the boundary. Thus the temperature distribution in R is modelled by the steady heat equation ∇ T =− in R, (5) where D is the diffusivity of the interior temperature T through the penguin huddle and R R Q is the (constant) source term representing the self-generation of heat by the penguins. In order to obtain parameter estimates for the interior problem, it is helpful to ﬁrst con- sider a circular penguin huddle of radius L, noting that T (0) = T . The Poisson equation R c (5) can then be solved exactly, giving Q Q 4(T − T ) c f T = T(r) = T − r , where = . (6) R c 4D D L R R The relation (6)isusedinSect. 2.1 to non-dimensionalise the problem for cases when the region R is non-circular. Finally, a condition on γ(t) is sought by specifying its outward normal velocity v .Fol- lowing [6], the idea that the boundary evolves owing to the heat loss of penguins on the boundary is used. Consequently, the problem has some similarity to the two-phase Stefan problem [23, 24, 29, 51] with the boundary temperature gradient driving its evolution, yet, unlike such problems, no phase change occurs. Waters et al. [6] move individual boundary penguins according to their heat loss to the exterior, whereas here both exterior and interior heat ﬂuxes are considered. First, the heat loss to the exterior wind is the ﬂux κ n ˆ ·∇T eval- uated on γ(t),where κ is the thermal conductivity of air. Second, the heat gain due to the warm huddle interior is −κ n ˆ ·∇T evaluated on γ(t),where κ is the thermal conductivity R R R of huddling penguins. Equating the normal velocity with the net heat ﬂux gives Av = κ n ˆ ·∇T − κ n ˆ ·∇T , on γ(t), (7) n R R where the constant term A is chosen to ensure a match in dimensions between both sides of (7), with the particular choice of A affecting the timescale for the boundary evolu- tion. For convenience in the non-dimensionalisation process upcoming in Sect. 2.1,let A = λρ c (T − T ),where ρ and c are the density and speciﬁc heat capacity of the f ∞ air, respectively, and λ is a non-dimensional constant to be chosen later. The right hand side of (7) is the jump in heat ﬂux at the boundary. This observation reinforces the analogy with the Stefan problem where typically the normal velocity of the boundary depends on the jump in the heat ﬂux [51]. For example, analogous boundary con- ditions arise in other natural free boundary problems, including the dissolution of ﬁnite 7 Page 6of17 S.J.Harris, N.R. McDonald objects, where only the gradient of the agent undergoing advection and diffusion exterior to the object is important [22, 24], or in porous media ﬂow about freeze pipes, in which the temperature gradient in both the exterior ﬂuid and interior ice regions are signiﬁcant [29]. However it is emphasised that, in this work, (7) is purely a modelling assumption. Rewriting (7) gives the following condition on γ ∂x λρ c (T − T ) n ˆ · = κ n ˆ ·∇T − κ n ˆ ·∇T + C(t). (8) f ∞ R R ∂t The additional term C (t) in (8) is needed to enforce conservation of penguins, or, equiva- lently, area conservation [6], namely that ∂x n ˆ · ds = v ds = 0. (9) ∂t γ γ 2.1 Non-dimensionalisation and Parameter Values In order to identify the parameters of the problem, the model (1)-(5), (8) is non- dimensionalised, with non-dimensional quantities temporarily labelled as starred variables. Let L characterise the length scale of region R i.e. x = Lx , and use the far-ﬁeld wind magnitude U as the characteristic velocity scale, u = U u . The following scalings are also used 1 T − T 1 U ∗ ∗ ∗ ∇ = ∇ ,T = ,t = t ,τ = . (10) L T − T τ L f ∞ Using (6), the non-dimensional form of the Poisson equation (5) can be written (now dropping the stars) as 4(T − T ) c f ∇ T =−α, α = . (11) T − T f ∞ The boundary condition (8) becomes ∂x λPe n ˆ · = n ˆ ·∇T − βn ˆ ·∇T + C(t), (12) ∂t where Pe = UL/D is the Péclet number, with D = κ /ρ c the thermal diffusivity of air, β = κ /κ is the ratio of thermal conductivities and C(t) is the dimensionless area conserving term, an explicit expression for which can be found by substituting (12)into(9) (see (13g) below). Note that the appearance of Pe on the left hand side of (12) affects the time scaling reported in the results. The real world dimensional time scale “t ” is therefore t = λPe (L/U )t . In what follows, the choice λ = 1 is made and justiﬁed a posteriori by comparing the time for the huddle to evolve with observations - see Sect. 5. There are three parameters: Pe, β and α. Following the justiﬁcation in [6], the wind around the penguin huddle is likely turbulent and hence the Péclet number can be approxi- mated by a turbulent Reynolds number, which is O(100) for Antarctic winds at this scale. Consistent with [6], the present work considers Péclet values between 1 and 100. The values of the remaining two parameters, α and β , are estimated as follows. The air temperature far upstream is taken to be T =−40 C[2, 3] while ambient temperatures in the huddle are above freezing [11]. It seems reasonable then to let T = 0 C. The centre of a huddle can f Penguin Huddling: A Continuum Model Page 7 of 17 7 ◦ ◦ exceed 20 C, with a maximum temperature of 37.5 C having been recorded [8]. Using, therefore, T = 30 Cgives α = 3. From table C.5 of [52], the thermal conductivity of air at −40 C is approximately κ ≈ −1 −1 −2 −1 0.022Wm K . The thermal conductance of penguins is 1.3Wm K from [1], however this is not the same as the thermal conductivity - it has different units for example. To ﬁnd the thermal conductivity, the conductance must be multiplied by the thickness of the −1 −1 penguin insulation layer which is d = 0.024m [53]. Therefore, κ = 0.0312Wm K and so β = 1.4182 ≈ 1.5 = O(1). This estimate for β depends on the wind ﬂow being either laminar or turbulent, with β decreasing as the ﬂow becomes more turbulent and the effective thermal conductivity of air increases [54]. In this work, a range of β values between β = 0 −1 and β = 1.5 are considered with β = 0.14 ≈ 0.1 = O(10 ) taken to be a representative value. To summarise, the full dimensionless system for the penguin problem is Pe v = n ˆ ·∇T − βn ˆ ·∇T + C(t), on γ, (13a) n R ∇ φ =0in , (13b) Pe u ·∇T =∇ T in , (13c) ∇ T =−α in R, (13d) ∂φ T = T = 1, =0on γ, (13e) ∂n T → 0,φ → x as r →∞, (13f) − (n ˆ ·∇T − βn ˆ ·∇T )ds C(t) = . (13g) ds Thetaskistosolve thesystem(13a)–(13g) and, in particular, ﬁnd the evolution of the γ(t) for some given initial huddle shape γ(0). This is done in Sect. 4 with the following choices: Pe = 1 → 100, β = 0 → 1.5and α = 3. 3 Numerical Method Following the work of [22, 29], the free boundary problem (13a)–(13g) can be written as a PG type equation [43] in terms of a conformal map z = f(ζ, t) from the exterior of the unit disk in the canonical ζ −plane to the exterior region in the physical z−plane [6, 26]. However, the interior R is not mapped using f(ζ, t), thus care must be taken with the term involving the interior temperature T . First, note the following which hold on the unit disk r =|ζ|= 1[20, 22] Re f ζf t ζ 1 ∂T σ(θ) v = , n ˆ ·∇T = Re[n∇T ]= = , (14) |f | |f | ∂r |f | ζ ζ ζ iθ where ∇= ∂ − i∂ , ζ = e , f = ∂f /∂ ζ and n = n + in = ζf /|f | is the complex x y ζ x y ζ ζ representation of the normal vector in the z−plane [34]. The function σ(θ) is found by solving the following integral equation [22, 26, 29] π = K Pe(cos θ − cos θ ) σ(θ )dθ , (15) 0 7 Page 8of17 S.J.Harris, N.R. McDonald where K(x) = e K (|x|) and K is the modiﬁed Bessel function. The function σ is the heat 0 0 ﬂux (due to the exterior effect alone) around the unit ζ −disk [22]. A numerical solution to (15) can be found at M equally spaced points on the upper half of the unit disk θ ∈ (0,π),see ∗ ∗ appendix A of [22]. To evaluate σ on the remainder of the disk, note that σ(θ ) = σ(π + θ ), where 0 ≤ θ <π,and σ(θ) → 1/π as θ → 0[22, 26]. Second, consider the contribution from the interior n ˆ ·∇T . This is found by direct com- putation in the z−plane. To begin, introduce the function u = T + αr /4; from (13d)and (13e), u satisﬁes ∇ u =0in R, (16) αr u = 1 + on γ. (17) This is a Laplace problem in the interior of a smooth, simply connected domain with Dirich- let boundary condition u = h(z) on γ . Efﬁcient methods developed by Trefethen and col- leagues are here used to solve the Laplace problem directly, see [47, 55] for further details. These are based on a rational approximation of u in the complex z−plane: there exists an analytic function F(z) such that u = Re[F(z)] which can be approximated by the series solution N N 1 2 F(z) = A + A (z − c ) + , (18) 0 k 0 z − p k=1 k=1 where A , A and B are constants, N is the power series truncation, c is a point in the 0 k k 1 0 interior R (in this work, the conformal centre of the polygon is used [22]) and p are some N poles in the exterior . The ﬁrst sum is known as the “Runge” (or smooth) part and the second as the “Newman” (or singular) part [49, 55]. Suitable poles p are found by the “adaptive Antoulas-Anderson” (AAA, pronounced ‘triple-A’) algorithm [49, 56] and the unknown coefﬁcients A , A and B by a least-squares 0 k k (LS) algorithm [46, 48, 57]: combining these into a AAA-LS method [50] can thus ﬁnd F(z) numerically. It follows [46]that ∇u =∇ Re[F(z)]= F (z), (19) where the primed notation indicates differentiation, in this case with respect to z. Therefore, ∇T = ∇[u − αr /4]= F (z) − z, (20) recalling that ∇= 2∂ and r = zz.Thisgives 1 ω(θ) n ˆ ·∇T = Re[n∇T ]= Re[ζf (F (z) − αz/2)]≡ , (21) R R ζ |f | |f | ζ ζ where the boundary points z in (21) are considered in the ζ −plane so that the right hand side of (21) is purely a function of θ . Drawing an analogue with σ , the function −βω is the heat ﬂux around the unit ζ −disk due to the interior effect alone. Third, the formula (13g) for the area conserving term C(t) can be transformed to the ζ −plane, noting that ds =|f |dθ . Substituting (14)and (21)into(13g)gives − σ − βω dθ −π C(t) = . (22) |f |dθ −π Penguin Huddling: A Continuum Model Page 9 of 17 7 Therefore, using (14), (21)and (22), the PG form of (13a)is Pe Re[f ζf ]= σ(θ) − βω(θ) + C(t)|f |, (23) t ζ ζ where σ(θ) satisﬁes (15)and ω(θ) = Re[ζf (F (z) − αz/2)], with F(z) given by (18). Note that σ(θ) − βω(θ) is the total heat ﬂux across |ζ|= 1; recall that the heat ﬂux across the huddle boundary γ is n ˆ ·∇T − βn ˆ ·∇T . It remains to ﬁnd the conformal map f , the general form of which is [24] −k z = f(ζ, t) = a (t )ζ + c (t )ζ , (24) −1 k k=0 where a (t ) is a real function in time and c (t ) = a (t ) + ib (t ) are complex functions −1 k k k in time. Note that a is the conformal radius and c is the conformal centre [22]. The −1 0 same c is used as the interior point in constructing the power series solution (18)for the Poisson problem in R. Numerically, the Laurent series (24) is truncated at N terms, giving n = 2(N + 1) + 1 = 2N + 3 unknown real functions in time: a (t ), a (t ), b (t ), k = −1 k k 0, 1,...,N . Likewise, n equally spaced points around the unit ζ −disk are chosen; the ﬁrst point is at θ = 0, the next N + 1 points are in the interval θ ∈ (0,π) and the ﬁnal N + 1 points in θ ∈ (π, 2π). Hence the choice M = N + 1in (0,π) is made when numerically solving the integral equation (15)for σ(θ). The PG equation (23) becomes a system of n coupled ordinary differential equations (ODEs) determining the time evolution of the n Laurent coefﬁcients (see [20, 21, 34]). In this work, the MATLAB routine ode15i was used to solve this ODE system. 4 Results In the results presented here, the evolution of the huddle boundary is plotted at equal time intervals in the range [0,t ],where t is the maximum time. The series truncation max max N = 128 is chosen, leading to a system of n = 259 ODEs to be solved. The initial start- ing shape of the penguin huddle is taken to be either a circle, a slanted ellipse, a triangle, an irregular pentagon or an hourglass. The initial Laurent coefﬁcients a (0), c (0) for the −1 k conformal map (24) of these shapes are taken from [24] (see their Figures 2a, 2b, 2d, 5a and 7a, respectively) and scaled such that each shape encloses the same area. Area conserva- tion of the huddle was monitored during each numerical experiment and found to hold with −5 a relative error of O(10 ). As a test, experiments were run which neglected the interior (i.e. β = 0) and the area conserving (C(t)) terms; the geometry of the obtained solutions matched quantitatively with results in [22] who considered an analogous problem in this class. First, consider the “exterior only” penguin problem when β = 0. This is a continuum ana- log of the discrete model developed in [6], where the penguin huddle evolves in response to the wind but the interior reorganisation of individual penguins is not represented. Figure 2 shows results for various Péclet values and initial starting shapes of the huddle. In general, these show good agreement to the experiments of [6]: the huddle propagates in the windward direction and evolves into a more streamlined shape over time. For large time (t t /2), max This scaling can be found analytically using the complex form of Green’s theorem [45] to ﬁnd the exact area from the conformal map z = f(ζ, 0). 7 Page 10 of 17 S.J. Harris, N.R. McDonald Fig. 2 Exterior only (β = 0) penguin problem. Shapes are plotted at equal time intervals up to a maximum time t .The x and y axes (not shown) are scaled at a 1-to-1 aspect ratio. Each description refers to the max initial penguin huddle shape, the Péclet number and the value of t : a) circle with Pe = 1, t = 10, b) max max circle with Pe = 10, t = 20, c) slanted ellipse with Pe = 10, t = 20, d) irregular pentagon with Pe max max = 100, t = 50 max the huddle approaches an invariant egg-shape which continues to propagate downwind with some constant velocity. To test this observation, a root mean-squared error (RMSE) measur- ing the deviation between consecutive plotted huddle shapes is computed for the examples −4 in Fig. 2 and found to approach a constant value with a relative error of O(10 ).Thisper- mits the reasonable conclusion that huddles do evolve into a steady shape and thus in all experiments presented, the maximum time t is chosen such that the penguin huddle has max reached a steady shape prior to or at time t . max Now, consider the full penguin problem where exterior, interior (β = 0) and area con- serving terms are all included - recall from Sect. 2.1 that β = 0.1 is used. Figures 3a-e present the evolution of the same initial starting shape, a circle, with Pe = 1, 5, 10, 50 and 100. The ﬁnal plot in each experiment is the steady shape and these are compared in Fig. 3f, where the shapes have been re-centered so that they share the same centre of mass. As the Péclet number is increased, the steady shape becomes less circular and more egg like. Similarly, the effect of changing the parameter β can be explored. Figures 4a-e shows the evolution of the circle with the values β = 0, 0.1, 0.5, 1 and 1.5, respectively, where Pe = 10 and t = 20 are taken for all plots. While a low value of β gives results that look max similar to the exterior only penguin problem, a high value of β causes the penguin huddle to retain a more circular shape. This is shown in Fig. 4f, where the steady shapes of Figs. 4a-e are compared. The dependence of the ﬁnal steady shape of the huddle on the values of the Péclet num- ber and β is explored further. A given steady shape can be quantiﬁed by its aspect ratio (AR): the ratio between its major and minor axes with AR = 1 being a circle and increasing AR > 1 becoming more egg-shaped. Figures 5a,b show the dependence of AR on Pe and β , respectively. Each marker on the diagram represents an experiment run: each curve contains 21 markers, all run to a maximum time of t = 80, which is sufﬁciently long to reach the max steady state huddle shape. The results from Fig. 5 agree with the conclusions drawn from Figs. 3 and 4.First,asPe increases, AR increases and thus a higher Péclet number results in a more egg-like steady shape, with a lower value of β giving higher values of AR for all time. This is also consistent with the parameter study in [6] - see their Fig. 3 - which shows that the huddle thickness decreases with an increasing Péclet number. As Pe → 0, AR→ 1, consistent with the limit of a circular steady shape when there is no wind, as expected. Second, as β increases, AR decreases, where lower values of Pe result in lower values of the AR for all time. Penguin Huddling: A Continuum Model Page 11 of 17 7 Fig. 3 The effect of different Péclet numbers on the evolution of the huddle is shown by plotting its shape at equal time intervals. The initial starting shape of the penguin huddle is a circle, with β = 0.1and a) Pe = 1, t = 10, b) Pe = 5, t = 15, c) Pe = 10, t = 20, d) Pe = 50, t = 40, e) Pe = 100, t = 50. f) max max max max max Comparison of the steady shapes of (a)-(e), the initial (circular) huddle shape is also plotted for comparison Fig. 4 Interior β effect: initial starting shape of a circle with Pe = 10, t = 20 and a) β = 0, b) β = 0.1, max c) β = 0.5, d) β = 1, e) β = 1.5. f) Comparison of the steady shapes of (a)-(e), the initial (circular) huddle shape is also plotted for comparison 7 Page 12 of 17 S.J. Harris, N.R. McDonald Fig. 5 The aspect ratio of the ﬁnal steady shape of the penguin huddle as a function of a) Péclet number, for the three given values β = 0, 0.1and 1, and b) β , for the three given values Pe = 10, 50 and 100. Each marker represents a numerical experiment Fig. 6 Initial shape effect. Shapes are plotted at equal time intervals with Pe = 10, β = 0.1, t = 35 and the max initial starting shape of a) a circle, b) a slanted ellipse, c) a triangle, d) an irregular pentagon, e) an hourglass. f) Comparison of the steady shapes of (a)-(e): all of the steady shapes are overlapped Figure 6 shows how huddles of different starting shapes evolve under the same Péclet number Pe = 10, the same interior effect β = 0.1 and up to the same maximum time t = max 35. Figures 6a-e show the evolution of a circle, a slanted ellipse, a triangle, an irregular pentagon and an hourglass, respectively, and Fig. 6f compares the steady shapes of these ﬁve experiments. The results show that, irrespective of the starting shape of the huddle and even those starting with signiﬁcant concave boundaries, the same steady shape is reached. The along-boundary variation of the total heat ﬂux across the huddle boundary of the steady shape of Fig. 6f is considered. By symmetry, only the upper half of the boundary need be considered; this is plotted in Fig. 7a. The arrows point in the normal direction of the Penguin Huddling: A Continuum Model Page 13 of 17 7 Fig. 7 a) Total heat ﬂux across the boundary of the steady shape of Fig. 6f. The arrows point in the normal direction of the given boundary point and their magnitude is proportional to the local heat ﬂux across the boundary. b) Canonical heat ﬂux in the upper half ζ −disk vs θ . c) Physical heat ﬂux in the upper half of the huddle boundary vs x given point on the boundary and their lengths represent the total physical heat ﬂux at that point. Recall from Sect. 3 that the total heat ﬂux in the physical domain across the huddle boundary γ is given by n ˆ ·∇T − βn ˆ ·∇T whereas the total heat ﬂux in the canonical domain across the unit ζ −disk is given by σ(θ) − βω(θ). The canonical heat ﬂux in the upper half ζ −disk vs θ is plotted in Fig. 7b and the physical heat ﬂux vs x is plotted in Fig. 7c. There is a much larger heat ﬂux on the windward side of the penguin huddle (the side of the penguin huddle directly facing the oncoming wind) compared to the leeward side. This is expected and is in agreement with the exterior temperature proﬁle presented in [6] - see their Fig. 1 - which shows that there is a steeper temperature gradient, and hence a higher total heat ﬂux, on the windward side. 5 Discussion In this continuum model, the evolution of the penguin huddle boundary γ(t), is affected by three factors: exterior cooling by the ambient wind, interior reorganisation of the penguins driven by their self-generation of heat, and area conservation of the huddle. The resulting two-dimensional model (13a)-(13g) can be rewritten as a PG-type equation (23)interms of the conformal map z = f(ζ, t) between the canonical (ζ ) and physical (z) planes. The map is approximated as a truncated Laurent series (24), whose unknown time-varying coefﬁcients are found by numerically solving a system of coupled ODEs. Results from Sect. 4 show that the huddle shape eventually reaches a steady shape. As the wind strength increases, the huddle becomes increasingly egg-shaped, minimising the portion of the huddle bound- ary directly facing the oncoming wind. For a stronger interior heating effect, the huddle 7 Page 14 of 17 S.J. Harris, N.R. McDonald approaches a more circle-like steady shape, to optimise interior heat regulation. For given wind and interior effects, any initial huddle shape eventually approaches the same steady shape. There is also a greater total heat ﬂux on the windward side of the huddle than on the leeward side. While dimensionless quantities have been used throughout the paper, it is useful to esti- mate the dimensional time scale associated with the huddle evolution. As noted in Sect. 2.1, ∗ ∗ the appropriate dimensional time scale is t = λPe(L/U )t where t is dimensionless time, and the choice λ = 1 has been made. A large huddle of, for example, 6000 penguins [1], is assumed and taking the radius of each penguin to be r ≈ 0.1m [3]gives L ≈ 10 m. Now −1 assume a sufﬁciently strong wind of U ≈ 20 ms for Pe = 100. Considering Fig. 3cwhere t = t = 50, the total real time for this experiment is t = 2500 seconds or ≈ 42 min- max r utes. Thus it is concluded that the penguin huddle motion described in this work occurs over roughly a one hour period, in agreement with previous modelling and ﬁeld observations [6, 8, 9], and justifying the choice λ = 1. In this paper, a continuum model has been used to ﬁnd the shape and propagation of the huddle boundary at the expense of tracking individual penguin movements. While there are pre-existing discrete models for penguin huddling (and other animal grouping dynam- ics) [4, 6, 7, 36], treating a penguin huddle as a continuum appears to be new. One of the advantages of the continuum model developed here is its computational efﬁciency, with a typical runtime of around 5-10 minutes on a standard laptop for a huddle to reach a steady shape. This can be sped up further, to a matter of seconds, by reducing the Laurent series truncation N , though care should be taken to ensure that the accuracy of the results is not compromised. When considering the “exterior only” (β = 0) penguin problem, our results matched well with those from [6], who did not include the interior effect of penguin reorganisation in their model. As discussed in Sect. 2.1, fully turbulent wind ﬂows have β ≈ 0 and hence are well modelled by considering exterior effects alone. The representative value β = 0.1, which is likely to be an overestimate, chosen in this work demonstrates that the interior heating of the colony only provides a weak effect in determining the steady shape of the huddle, see Fig. 4f. Larger values of β resulting in more circular huddle shapes have been included here for completeness. Further, the results in this work match well with real world observations: penguin huddles are seen to march downwind over time and observed huddle shapes resemble the egg-like steady shapes found here - see [1, 4, 10, 12, 13]. Extensions to the penguin problem are of interest. A time-varying wind could be consid- ered, where the variations are on the order of an hour; when the wind direction changes, the huddle will seek to propagate back to the optimal (egg) steady shape in the new windward direction. Stochastic terms could also be incorporated into the model, as in [6]. This includes random variations in the direction and magnitude of the wind and stochasticity in the interior and exterior heat ﬂuxes. The numerical method would need to be adapted to handle these stochastic terms. The evolution of any arbitrary starting huddle shape could also be simulated. In this work, the initial huddle boundary was taken such that its conformal map was readily represented in the form of a truncated Laurent series (24). For more arbitrary starting shapes, the hud- dle boundary can be discretised as a polygon of m vertices, then its conformal map to the ζ −plane found numerically at each time step as the penguin huddle evolves. Pre-exisiting numerical schemes for ﬁnding the conformal map between a polygon and the unit ζ −disk include the Schwarz-Christoffel toolbox [58] and a variation of the AAA-LS algorithm [50] from Sect. 3. It would also be possible to consider multi-connected huddles such as those with penguin-free holes or the interaction, and possible merger, of initially disjoint huddles. Penguin Huddling: A Continuum Model Page 15 of 17 7 While this would require a signiﬁcant modiﬁcation of the numerical method used here, it is worth remarking that AAA-LS type methods are able to handle such geometries [47, 49, 50]. Finally, consider the boundary condition (13a), and its associated PG equation (23). The right-hand side is constructed by exterior (n ˆ ·∇T ), interior (βn ˆ ·∇T ) and area conserva- tion (C(t)) effects. By inclusion and exclusion of these three terms, (13a) encapsulates six possible free boundary problems: 1. The dissolution problem: exterior effect only. A permeable object is placed in a uniform stream of some dissolving agent - see [22, 59]. 2. The Poisson growth problem: interior effect only. Applicable to squeeze ﬂow and the evaporation of thin liquid ﬁlms - see [45, 60, 61]. 3. The two-phase melting/freezing problem: exterior and interior effects. Porous media ﬂow about freeze pipes where interior and exterior temperature gradients are in effect - see [24, 29]. 4. The “exterior only” penguin problem: exterior and area conservation effects. Huddle problems where interior reorganisation is not represented, such as if individuals in the centre are unable to move or there exists a “pecking order” in the population - see [6]. 5. The “bat huddle” problem: interior and area conservation effects. Huddle problems with no wind effect, for example bats huddling in caves sheltered from exterior winds. To the authors’ knowledge, there is no mathematical work on this problem, but for biological context see [39, 40]. 6. The “full” penguin problem: exterior, interior and area conservation effects. This is the full huddle continuum problem covered in this work. Therefore, the numerical method developed in this work can be used to model and simulate additional free boundary problems, such as those from the ﬁelds of ﬂuid mechanics and mathematical biology. Acknowledgements The authors thank Carole Hall, Stony Brook University, for their inspiring talk on Adélie penguins and the useful discussion thereafter. The authors also thank the anonymous reviewers for their useful comments and suggestions. Author Contribution Both authors contributed equally to this work. Funding Samuel J. Harris was supported by a UK Engineering and Physical Sciences Research Council PhD studentship, grant numbers EP/N509577/1 and EP/T517793/1. Data Availability All data generated or analysed during this study are included in this published article and its supplementary information ﬁles. Code Availability All relevant code ﬁles are publicly available at the following GitHub link: https://github. com/Sam-J-Harris/Penguin-huddling-a-continuum-model.git. Declarations Competing Interests Not applicable. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 7 Page 16 of 17 S.J. Harris, N.R. McDonald References 1. Le Maho, Y.: The emperor penguin: a strategy to live and breed in the cold: morphology, physiology, ecology, and behavior distinguish the polar emperor penguin from other penguin species, particularly from its close relative, the king penguin. Am. Sci. 65(6), 680–693 (1977) 2. McCafferty, D.J., Gilbert, C., Thierry, A.-M., Currie, J., Le Maho, Y., Ancel, A.: Emperor penguin body surfaces cool below air temperature. Biol. Lett. 9(3), 20121192 (2013) 3. Williams, C.L., Hagelin, J.C., Kooyman, G.L.: Hidden keys to survival: the type, density, pattern and functional role of emperor penguin body feathers. Proc. R. Soc. Lond. B 282(1817), 20152033 (2015) 4. Gerum, R.C., Fabry, B., Metzner, C., Beaulieu, M., Ancel, A., Zitterbart, D.P.: The origin of traveling waves in an emperor penguin huddle. New J. Phys. 15(12), 125022 (2013) 5. Kooyman, G.L., Gentry, R.L., Bergman, W.P., Hammel, H.T.: Heat loss in penguins during immersion and compression. Comp. Biochem. Physiol. 54(1), 75–80 (1976) 6. Waters, A., Blanchette, F., Kim, A.D.: Modeling huddling penguins. PLoS ONE 7(11), 50277 (2012) 7. Gerum, R., Richter, S., Fabry, B., Le Bohec, C., Bonadonna, F., Nesterova, A., Zitterbart, D.P.: Structural organisation and dynamics in king penguin colonies. J. Phys. D, Appl. Phys. 51(16), 164004 (2018) 8. Gilbert, C., Robertson, G., Le Maho, Y., Naito, Y., Ancel, A.: Huddling behavior in emperor penguins: dynamics of huddling. Physiol. Behav. 88(4–5), 479–488 (2006) 9. Kirkwood, R., Robertson, G.: The occurrence and purpose of huddling by emperor penguins during foraging trips. Emu 99(1), 40–45 (1999) 10. Ancel, A., Gilbert, C., Poulin, N., Beaulieu, M., Thierry, B.: New insights into the huddling dynamics of emperor penguins. Anim. Behav. 110, 91–98 (2015) 11. Zitterbart, D.P., Wienecke, B., Butler, J.P., Fabry, B.: Coordinated movements prevent jamming in an emperor penguin huddle. PLoS ONE 6(6), 20260 (2011) 12. Mina, T., Min, B.-C.: Penguin huddling inspired distributed boundary movement for group survival in multi-robot systems using Gaussian processes. In: 2018 IEEE International Conference on Robotics and Biomimetics, pp. 2177–2183. IEEE (2018) 13. Richter, S., Gerum, R., Winterl, A., Houstin, A., Seifert, M., Peschel, J., Fabry, B., Le Bohec, C., Zit- terbart, D.P.: Phase transitions in huddling emperor penguins. J. Phys. D, Appl. Phys. 51(21), 214002 (2018) 14. Gu, W., Christian, J.K., Woodson, C.B.: A novel coupled ﬂuid-behavior model for simulating dynamic huddle formation. PLoS ONE 13(8), 0203231 (2018) 15. Saffman, P.G., Taylor, G.I.: The penetration of a ﬂuid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. A 245(1242), 312–329 (1958) 16. Howison, S.: Fingering in Hele-Shaw cells. J. Fluid Mech. 167, 439–453 (1986) 17. Paterson, L.: Radial ﬁngering in a Hele Shaw cell. J. Fluid Mech. 113, 513–529 (1981) 18. Taylor, G., Saffman, P.G.: A note on the motion of bubbles in a Hele-Shaw cell and porous medium. Q. J. Mech. Appl. Math. 12(3), 265–279 (1959) 19. Entov, V.M., Etingof, P.I.: Bubble contraction in Hele-Shaw cells. Q. J. Mech. Appl. Math. 44(4), 507–535 (1991) 20. Dallaston, M.C., McCue, S.W.: An accurate numerical scheme for the contraction of a bubble in a Hele– Shaw cell. ANZIAM J. 54, 309–326 (2013) 21. Dallaston, M.C., McCue, S.W.: A curve shortening ﬂow rule for closed embedded plane curves with a prescribed rate of change in enclosed area. Proc. R. Soc. A 472(2185), 20150629 (2016) 22. Ladd, A.J.C., Yu, L., Szymczak, P.: Dissolution of a cylindrical disk in Hele-Shaw ﬂow: a conformal- mapping approach. J. Fluid Mech. 903, 46 (2020) 23. Cummings, L.M., Hohlov, Y.E., Howison, S.D., Kornev, K.: Two-dimensional solidiﬁcation and melting in potential ﬂows. J. Fluid Mech. 378, 1–18 (1999) 24. Rycroft, C.H., Bazant, M.Z.: Asymmetric collapse by dissolution or melting in a uniform ﬂow. Proc. R. Soc. A 472(2185), 20150531 (2016) 25. Mullins, W.W., Sekerka, R.F.: Stability of a planar interface during solidiﬁcation of a dilute binary alloy. J. Appl. Phys. 35(2), 444–451 (1964) 26. Choi, J., Margetis, D., Squires, T.M., Bazant, M.Z.: Steady advection–diffusion around ﬁnite absorbers in two-dimensional potential ﬂows. J. Fluid Mech. 536, 155–184 (2005) 27. Tsai, V.C., Wettlaufer, J.S.: Star patterns on lake ice. Phys. Rev. E 75(6), 066105 (2007) 28. Grodzki, P., Szymczak, P.: Reactive-inﬁltration instability in radial geometry: from dissolution ﬁngers to star patterns. Phys. Rev. E 100(3), 033108 (2019) 29. Goldstein, M.E., Reid, R.L.: Effect of ﬂuid ﬂow on freezing and thawing of saturated porous media. Proc. R. Soc. A 364(1716), 45–73 (1978) 30. Langer, J.S.: Instabilities and pattern formation in crystal growth. Rev. Mod. Phys. 52(1), 1 (1980) Penguin Huddling: A Continuum Model Page 17 of 17 7 31. Brower, R.C., Kessler, D.A., Koplik, J., Levine, H.: Geometrical models of interface evolution. Phys. Rev. A 29(3), 1335 (1984) 32. Sethian, J.A.: Curvature and the evolution of fronts. Commun. Math. Phys. 101(4), 487–499 (1985) 33. Hilton, J.E., Miller, C., Sharples, J.J., Sullivan, A.L.: Curvature effects in the dynamic propagation of wildﬁres. Int. J. Wildland Fire 25(12), 1238–1251 (2016) 34. Harris, S.J., McDonald, N.R.: Fingering instability in wildﬁre fronts. J. Fluid Mech. 943, 34 (2022) 35. Burger, M., Haškovec, J., Wolfram, M.-T.: Individual based and mean-ﬁeld modeling of direct aggrega- tion. Physica D 260, 145–158 (2013) 36. Bernardi, S., Scianna, M.: An agent-based approach for modelling collective dynamics in animal groups distinguishing individual speed and orientation. Philos. Trans. R. Soc. Lond. B 375(1807), 20190383 (2020) 37. Katz, Y., Tunstrøm, K., Ioannou, C.C., Huepe, C., Couzin, I.D.: Inferring the structure and dynamics of interactions in schooling ﬁsh. Proc. Natl. Acad. Sci. 108(46), 18720–18725 (2011) 38. Bhattacharya, K., Vicsek, T.: Collective decision making in cohesive ﬂocks. New J. Phys. 12(9), 093019 (2010) 39. Herreid, C.F.: Temperature regulation of Mexican free-tailed bats in cave habitats. J. Mammal. 44(4), 560–573 (1963) 40. Ryan, C.C., Burns, L.E., Broders, H.G.: Changes in underground roosting patterns to optimize energy conservation in hibernating bats. Can. J. Zool. 97(11), 1064–1070 (2019) 41. Nave, G.K. Jr, Mitchell, N.T., Chan Dick, J.A., Schuessler, T., Lagarrigue, J.A., Peleg, O.: Attraction, dynamics, and phase transitions in ﬁre ant tower-building. Front. Robot. AI 7, 25 (2020) 42. Ko, T.-Y.H., Yu, H.D.L.: Fire ant rafts elongate under ﬂuid ﬂows. Bioinspir. Biomim. 17(4), 045007 (2022) 43. Bazant, M.Z., Crowdy, D.: Conformal mapping methods for interfacial dynamics. In: Handbook of Ma- terials Modeling, pp. 1417–1451. Springer, Dordrecht (2005) 44. Gustafsson, B., Vasil’ev, A.: Conformal and Potential Analysis in Hele-Shaw Cells. Springer, Berlin (2006) 45. McDonald, R., Mineev-Weinstein, M.: Poisson growth. Anal. Math. Phys. 5(2), 193–205 (2015) 46. Trefethen, L.N.: Series solution of Laplace problems. ANZIAM J. 60(1), 1–26 (2018) 47. Trefethen, L.N.: Numerical conformal mapping with rational functions. Comput. Methods Funct. Theory 20(3), 369–387 (2020) 48. Baddoo, P.J.: Lightning solvers for potential ﬂows. Fluids 5(4), 227 (2020) 49. Costa, S.: Solving Laplace problems with the AAA algorithm (2020). arXiv:2001.09439. ArXiv preprint 50. Costa, S., Trefethen, L.N.: AAA-least squares rational approximation and solution of Laplace problems (2021). arXiv:2107.01574. ArXiv preprint 51. Gupta, S.C.: The Classical Stefan Problem: Basic Concepts, Modelling and Analysis with Quasi- Analytical Solutions and Methods, vol. 45. Elsevier, Amsterdam (2017) 52. Basu, P.: Biomass Gasiﬁcation, Pyrolysis and Torrefaction: Practical Design and Theory. Academic Press, Boston (2018) 53. Dawson, C., Vincent, J.F.V., Jeronimidis, G., Rice, G., Forshaw, P.: Heat transfer through penguin feath- ers. J. Theor. Biol. 199(3), 291–295 (1999) 54. Kittel, C., McEuen, P., McEuen, P.: Introduction to Solid State Physics, vol. 8. Wiley, New York (1996) 55. Gopal, A., Trefethen, L.N.: Solving Laplace problems with corner singularities via rational functions. SIAM J. Numer. Anal. 57(5), 2074–2094 (2019) 56. Nakatsukasa, Y., Sète, O., Trefethen, L.N.: The AAA algorithm for rational approximation. SIAM J. Sci. Comput. 40(3), 1494–1522 (2018) 57. Brubeck, P.D., Nakatsukasa, Y., Trefethen, L.N.: Vandermonde with Arnoldi. SIAM Rev. 63(2), 405–415 (2021) 58. Driscoll, T.A.: Algorithm 756: a MATLAB toolbox for Schwarz-Christoffel mapping. ACM Trans. Math. Softw. 22(2), 168–186 (1996) 59. Dutka, F., Starchenko, V., Osselin, F., Magni, S., Szymczak, P., Ladd, A.J.: Time-dependent shapes of a dissolving mineral grain: comparisons of simulations with microﬂuidic experiments. Chem. Geol. 540, 119459 (2020) 60. Agam, O.: Viscous ﬁngering in volatile thin ﬁlms. Phys. Rev. E 79(2), 021603 (2009) 61. Crowdy, D., Kang, H.: Squeeze ﬂow of multiply-connected ﬂuid domains in a Hele-Shaw cell. J. Non- linear Sci. 11(4), 279–304 (2001) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afﬁliations.
Acta Applicandae Mathematicae – Springer Journals
Published: Jun 1, 2023
Keywords: Free boundary problems; Conformal mapping; AAA-least squares algorithm; Fluid mechanics; Collective behaviour; Huddling
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.