A second example of the RFID system allows for tests of spatial cognition (spatial learning and memory) in Mountain Chickadees (Poecile gambeli) in the Sierra Nevada Mountains. These birds live year-round in mountainous habitat that can receive up to six meters of snow depth in the winter. To survive these conditions the birds must store or cache food (usually pine seeds) and then recover them days or weeks later when other food sources are unavailable in winter. Hence, a large capacity for spatial cognitive ability is key to the survival of these birds, and individual variation in spatial learning and memory ability is of great interest in understanding the ecology and evolution of these animals as well as food-caching species in general.

Researchers at the University of Nevada (VP, CB, DK, AP, and BS) have devised and implemented a means of evaluating spatial learning and memory in a population of Mountain Chickadees using arrays of RFID-enabled bird feeders that can selectively feed specific individuals. The feeders are equipped with a motorized door, controlled by the RFID circuit board, that can be lowered or raised to allow or deny access to food. An antenna is embedded in a perch positioned in front of the door such that the feeder can recognize individuals on the perch and operate the door accordingly to provide or withhold food. Movement of the door is accomplished with a rack-and-pinon gear mechanism powered by a small gear motor. A limit switch is used in conjunction with the door to indicate to the circuit board when the door had reached a fully open or fully closed position. The motor is controlled by a customized accessory circuit board that features a TB6612FNG motor controller (details on OSF project page).

The learning and memory tests employed arrays of eight RFID-enabled bird feeders that were mounted together on a square frame with two feeders on each side facing outward (see **Figure 2B** and Pitera et al., 2018). One or two of these arrays (depending on the study) were situated at high and low elevations sites. The arrays were suspended on cables stretched between trees so the feeders could be raised above the reach of bears and rodents. Prior to testing, the feeders were all configured in the open position. That is, the doors were all open such that all birds have access to clearly visible food (sunflower seeds). RFID data logging was enabled at this time such that the researchers could determine which tagged birds were visiting the feeders, but the RFID reads did not affect access to food.

After this initial acclimation period of about a week, the feeders were switched to "feed-all mode," wherein the door remained closed until any bird with a tag is detected on the perch. In "feed-all" mode, any bird with a tag will cause the feeder door to open. This training mode allowed birds to become accustomed to the movement of the door.

After another initialization period in feed-all mode, the feeders were reconfigured into "target mode" such that each bird will have access to food at only one of the eight feeders. Spatial learning and memory were assessed based on how many nonrewarding feeders a bird visited prior to visiting the assigned, rewarding feeder during each trial. A trial starts when a bird visits any feeder in the array and ends with the visit to the rewarding feeder. Birds that quickly learn and remember the location of the rewarding feeder (as evinced by progressively fewer visits to non-rewarding feeders) are deemed as having superior spatial learning and memory ability. In addition to testing spatial cognition, the system can test reversal spatial learning and memory performance by switching the rewarding feeder for each bird after the completion of the spatial learning and memory task (which usually took 4 days). The number of errors (e.g., number of non-rewarding feeders visited prior to visiting the rewarding feeder) during a 4-days reversal trial provides a measure of spatial learning and memory flexibility as a bird needs to stop visiting the feeder that provided food in the previous spatial learning and memory task and learn the location of a new rewarding feeder.

This system was applied to Mountain Chickadees living at two different elevations. The birds at high elevation (∼2,400 m) face more severe and longer winter conditions and rely more heavily on caching and recovering food for overwinter survival than do birds at lower elevation (∼1,900 m). The RFID system provided a means of testing several key hypotheses relating spatial cognitive ability to environmental conditions and fitness. First, the spatial cognitive tests revealed that the high elevation birds had better spatial learning and memory ability than did the low elevation birds (Croston et al., 2016). Second, the system revealed that individual variation in spatial learning and memory performance is associated with differences in survival in first-year, juvenile birds during their first winter, showing that spatial cognition is affected by natural selection at high elevations (Sonnenberg et al., 2019). Moreover, a series of spatial learning and memory reversal tests suggested a potential tradeoff between cognitive flexibility and spatial learning/memory (Croston et al., 2017; Tello-Ramos et al., 2018). RFID data collected in this system were also used to show that daily foraging routines in chickadees differ between elevations and among seasons. Moreover, it was apparent that spatial learning and memory performance was associated with daily foraging routines. In particular, chickadees with better spatial cognition had daily foraging routines that resembled those in milder environment and seasons, likely due to greater predictability of foraging success for these individuals (Pitera et al., 2018). Most recently the system has revealed that females allocate more reproductive effort when mated with males with better spatial learning and memory performance (Branch et al., 2019). Finally, ongoing work in this system involves analyzing chickadee social networks using RFID data to address the associations between social and cognitive phenotypes as well as the potential role of cognitive phenotype in structuring these social networks.

### Implementation 3: A Nest Box Platform for Experimental Manipulation

Several studies have already employed RFID to generate a detailed activity log for birds that use artificial nest boxes (e.g., Johnson et al., 2013; Stanton et al., 2016; Zarybnicka et al., 2016; Schuett et al., 2017; Chien and Chen, 2018; Firth et al., 2018). With this third application example we describe a project that takes the next logical step forward—using RFID to orchestrate experimental treatments that can be applied individually to breeding adults and offspring. This project involved a nest box designed for Eastern Bluebirds (Sialia sialis) that could manipulate environmental noise and/or nestbox temperature in accordance with which birds were inside the box (**Figure 2C**). The project has progressed through several versions of the RFIDenabled nestbox, but not all features have been tested in the field (notably the temperature manipulation). Hence, the functions described here should be regarded as features that are available for a technology-driven nestbox study.

To help ensure that treatments are applied accurately to the targeted individuals, the nestbox was configured to make use of two antennas. One antenna was mounted to encircle the entrance to the nest box, and it would capture birds as they entered or exited. The second antenna was positioned inside the nest box, and it could verify the presence of a bird therein. The tag-polling strategy was conceived such that the entrance antenna was used for most of the monitoring, and polling by the internal antenna was only done periodically or when the entrance antenna had received a tag. This arrangement was necessary because birds will often rest on the nestbox entrance but not go in. Or they may fly up to the entrance but quickly fly away. The dual antenna system provided assurance that a bird had entered the nestbox.

Following the general example of Lendvai et al. (2015b), artificial noise was incorporated into the system by means of a generic serial MP3 music player module and a 1 W speaker (catalex.taobao.com). These items were purchased as a single kit for about \$4. The MP3 player accesses and plays audio tracks from a mini-SD card. The microcontroller on the RFID circuit board communicated with the MP3 player via a simple serial interface. Hence, the RFID reader could control the timing, duration, and volume of playbacks. We used a transistor to control the power supply to the MP3 player, such that it could be powered down when it was not being used (see OSF project page for details). Similarly, we have configured the RFID circuit board to control the power supply to an electronic heating pad (WireKinetics Co, Ltd, Taipei, Taiwan). As with the speaker system we used a simple transistor circuit to provide current to the heating pad directly from the power supply (i.e., battery).

As part of a pilot study, a single noise emitting system was deployed on an active bluebird nest in Kent County Michigan. The system was programmed to emit noise from 05:30 to 11:00 each day when both parents were absent from the nest box. The system was successful in administering an effective treatment throughout the majority of a nesting effort. The sample size is, of course, too small to reach any conclusions, and based on this initial effort, we plan to continue this experiment in future breeding seasons. The heating pad has not yet been tested, but we plan to do so in the spring of 2019, and the results will be posted to our OSF project page.

### DISCUSSION

The example studies described above do not by any means exhaust the possibilities offered by our low-cost RFID system. Given the wide array of hardware and software tools available as part of the Arduino electronics platform there are many other applications with respect to data logging and experimental manipulations that could be brought to bear in the context of RFID-based biologging. A few potential and/or recently executed applications are as follows:

As new applications come online, we intend to incorporate them as modules into our OSF project page.

There is a clear taxonomic bias in the examples described in this paper; they are all applied to avian study systems. However, this bias is due to the common interests of the research teams and not limitations in the applicability of RFID systems to other taxa. PIT tags have been used in a wide variety of vertebrates and even on a few insects (Sumner et al., 2007; Bonter and Bridge, 2011; Dogan et al., 2016; Whitham and Miller, 2016). Nevertheless, their applicability to bird research should not be surprising majority of bird species require small (<1 g) tracking devices.

A major advantage to automated systems like the ones described here is the potential for continuous data recording. Traditional monitoring methods often involve sampling the activity schedules of animals—for example, a researcher may visually monitor a burrow entrance for 3 h every other day, or one may video record activity at a feeding station for 5 h every day. Although these methods may be adequate for finding averages associated with events that happen regularly, they may be inadequate for detecting rare or even infrequent events (like brood parasitism or fledging). Automated systems obviate the need for sampling such that you can capture the entirety of a breeding attempt or activity period (see Lendvai et al., 2015a).

As is the case with all research efforts that involve attaching devices to animals, researchers must make every effort to minimize the adverse effects of equipping animals with PIT tags. There are numerous studies of the effects of PIT tags on survival and behavior, with most of this research focused on fish (Keck, 1994; Low et al., 2005; Nicolaus et al., 2008; Burdick, 2011; Thorstad et al., 2013; Guimaraes et al., 2014; Ratnayake et al., 2014; Moser et al., 2017; Schlicht and Kempenaers, 2018). The vast majority of these studies conclude that the effect of implanted PIT tags is negligible. The Wood Duck research described above entailed injecting tags in both females and ducklings intrascapularly using a PIT-tag syringe, and in 6 years, no deleterious effects have been detected on females or ducklings in the wild. Related projects have involved tagging >500 ducklings in captivity and raising >200 to the fledging stage with no adverse effects on survival or health observed in the course of regular monitoring by staff veterinarians. Tag loss has been <1% (Eadie et al., in preparation).

We do not know of a systematic study that has investigated the effects of externally mounted PIT tags. In particular, external mounting would largely apply to bird species, wherein PIT tags are incorporated into leg bands (see **Figure 1**). Our own experience has indicated no problems other than those typically associated with leg bands. More specifically, the studies of Mountain Chickadees described above have involved tagging over 1,000 individuals with RFID leg bands. Thus far, there has been one instance where the leg band appeared to cause an injury. However, there are unpublished reports of apparent injuries and mortality in association with RFID leg bands (Curry, unpublished data; Morand-Ferron et al., personal communication). We advocate the use of compact RFID leg bands like those supplied by Eccel technology, which minimize the size of the tag (see **Figure 1**), and we recommend caution in selecting the appropriate band size for the species (or individual) being tagged.

It is important to note that our device is not the only RFID reader available for field biologists. Among the alternatives for low-frequency (120–150 MHz) readers, the lowest-cost options include several RFID reader modules that can simply read tag data and communicate with a computer or microcontroller. Examples include the ID-12 and ID-20 readers (ID Innovations, Canning Vale, WA, Australia), the Paralax RFID module (Parallax, Inc, Rocklin, CA, USA), the MIKROE-262 and MIKROE-1434 modules (Mikroelektronika D.O.O., Belgrade, Serbia), and the RFIDREAD-RW Module (Priority 1 Design, Melbourne, Australia). These sorts of modules range in cost from about \$10–50 (USD). There are a few commercially available RFID readers that can log data in a manner similar to our ETAG reader. For example, Priority 1 Design offers the RFIDLOG circuitboard for ∼\$50 (USD). A more robust, fully enclosed RFID reader/datalogger is available from Eccel Technology, Ltd (Leicester, UK), for ∼\$400 (USD). For applications that require read ranges >2–3 cm, there are high power RFID readers capable of read distances of 50 cm or more. Unfortunately, readers that offer this increased read range will typically cost considerably more than the low power systems. Examples of high-power systems include Biomark readers, such as the IS1001, which costs ∼\$1,500 (USD), and systems from Oregon RFID (Portland, OR, USA) with a minimum cost of just over \$2,000 (USD). Although there are many RFID readers available, there are few if any standalone systems that can be programmed by the user to carry out complex protocols. Our ETAG reader offers this key feature.

To facilitate widespread collaboration and sharing of information, we have established a project within the Open Science Framework at https://osf.io/9j7ax/. This information repository contains technical information relating to the RFID system and the research applications described in this paper, and it is open to the public to view and download materials. It also allows for ongoing updates and expansion such that we can continue to add new applications and contributors into the foreseeable future. New developments on the horizon include wireless networks for real-time data delivery, and online software tools for designing antennas that will culminate in an online simulator for testing virtual antennas prior to making a physical prototype. We are also developing web-based data portal and archive for uploading, storing and managing RFID data. We hope that our RFID system and the forthcoming tools will prove to be useful resources for the biologging community, and we invite our readers to take part in our efforts to expand the utility of RFID technology for animal tracking and monitoring.

### ETHICS STATEMENT

All animal research was conducted according to accepted ethical standards and with approval from the authors' respective Institutional Animal Care and Use Committees (IACUC). Specifically, the research described in this manuscript followed IACUC-approved protocols R16-010, 00603, and 20971 for the University of Oklahoma, the University of Nevada Reno, and the University of California Davis, respectively. All field efforts sought to minimize negative impacts on birds by working with animals only during favorable conditions and by following best practices for animal handling and tagging.

### AUTHOR CONTRIBUTIONS

EB participated in designing and implementing the RFID systems and drafted the main text of the article. JW and AM helped design the RFID circuit board, wrote key elements of the firmware, and served as a grant PI. MP and DP designed and implemented the experimental bluebird nestbox and contributed to writing this section of the text. CH designed and implemented the solar charging system used on the bluebird boxes. CC served as a grant PI and along with TP has built data infrastructure for the RFID system and managed the OSF project page. JE, TS, AO, and BL carried out the field study of Wood Ducks, and JE and BL helped draft this section of text. CB, AP, DK, BS, and VP conceived and executed the spatial memory testing. VP and AP helped to draft this section of the text. JR served as the overall project manager and grant PI for the project, and she provided editorial assistance on the manuscript.

#### FUNDING

Funding for the technology development came from the National Science Foundation: NSF IDBR 1556313 to JR and EB; NSF IDBR 1556316 to JW; NSF ABI 1458402 to CC, JR, and EB. The NSF also provided funding for two of the featured RFID implementations: IOS 1351295 and NSF IOS 1856181 to VP; NSF IOS 1355208 to JE and BL. Development and testing of the bluebird nestbox was supported by a George Miksch Suttton Scholarship, a grant from the North American Bluebird Society, and an OU Biological Station Award to MP.

#### ACKNOWLEDGMENTS

We wish to thank all of our colleagues who tested prototypes and provided feedback. We also thank the many field

#### REFERENCES

assistants and student volunteers who were involved in the field research described in this paper. Finally, thanks also to Elena Bridge for assembling hundreds of RFID circuit boards.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Bridge, Wilhelm, Pandit, Moreno, Curry, Pearson, Proppe, Holwerda, Eadie, Stair, Olson, Lyon, Branch, Pitera, Kozlovsky, Sonnenberg, Pravosudov and Ruyle. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Dynamic-Parameter Movement Models Reveal Drivers of Migratory Pace in a Soaring Bird

Joseph M. Eisaguirre1,2 \*, Marie Auger-Méthé3,4, Christopher P. Barger 5 , Stephen B. Lewis 6 , Travis L. Booms 5 and Greg A. Breed1,7

1 Department of Biology and Wildlife, University of Alaska Fairbanks, Fairbanks, AK, United States, 2 Department of Mathematics and Statistics, University of Alaska Fairbanks, Fairbanks, AK, United States, 3 Department of Statistics, University of British Columbia, Vancouver, BC, Canada, 4 Institute for the Oceans and Fisheries, University of British Columbia, Vancouver, BC, Canada, 5 Alaska Department of Fish and Game, Fairbanks, AK, United States, 6 United States Fish and Wildlife Service, Juneau, AK, United States, 7 Institute of Arctic Biology, University of Alaska Fairbanks, Fairbanks, AK, United States

#### Edited by:

Frants Havmand Jensen, Woods Hole Oceanographic Institution, United States

#### Reviewed by:

Anne K. Scharf, Max Planck Institute of Ornithology, Germany Peter E. Smouse, Rutgers University, The State University of New Jersey, United States

> \*Correspondence: Joseph M. Eisaguirre<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 11 March 2019 Accepted: 07 August 2019 Published: 27 August 2019

#### Citation:

Eisaguirre JM, Auger-Méthé M, Barger CP, Lewis SB, Booms TL and Breed GA (2019) Dynamic-Parameter Movement Models Reveal Drivers of Migratory Pace in a Soaring Bird. Front. Ecol. Evol. 7:317. doi: 10.3389/fevo.2019.00317 Long distance migration can increase lifetime fitness, but can be costly, incurring increased energetic expenses and higher mortality risks. Stopover and other en route behaviors allow animals to rest and replenish energy stores and avoid or mitigate other hazards during migration. Some animals, such as soaring birds, can subsidize the energetic costs of migration by extracting energy from flowing air. However, it is unclear how these energy sources affect or interact with behavioral processes and stopover in long-distance soaring migrants. To understand these behaviors and the effects of processes that might enhance use of flight subsidies, we developed a flexible mechanistic model to predict how flight subsidies drive migrant behavior and movement processes. The novel modeling framework incorporated time-varying parameters informed by environmental covariates to characterize a continuous range of behaviors during migration. This model framework was fit to GPS satellite telemetry data collected from a large soaring and opportunist foraging bird, the golden eagle (Aquila chrysaetos), during migration in western North America. Fitted dynamic model parameters revealed a clear circadian rhythm in eagle movement and behavior, which was directly related to thermal uplift. Behavioral budgets were complex, however, with evidence for a joint migrating/foraging behavior, resembling a slower paced fly-andforage migration, which could facilitate efficient refueling while still ensuring migration progress. In previous work, ecological and foraging conditions are usually considered to be the key aspects of stopover location quality, but taxa, such as the golden eagle, that can tap energy sources from moving fluids to drive migratory locomotion may pace migration based on both foraging opportunities and available flight subsidies.

Keywords: Bayesian, correlated random walk, golden eagle, movement ecology, soaring flight

## INTRODUCTION

Long-distance migration can relax competition and permit use of seasonally available resources, helping many animals maximize lifetime fitness (Newton, 2008; Avgar et al., 2014). Those benefits, however, come at substantial costs, including greater vulnerability to predators, uncertain conditions, mechanical wear, elevated energy expenditure, and time (Alerstam and Hedenström, 1998; Clark and Butler, 1999; Hedenström, 2008; Newton, 2008; Avgar et al., 2014). As many migrant species cannot store sufficient energy for nonstop, long-distance migration, stopover evolved as a behavior for strategically resting and refueling en route (Gill, 2007).

Migrant species are adapted for utilizing either soaring or flapping flight, and the different flight modes translate into stopover strategy (Hedenström, 1993; Gill, 2007). Generally, soaring flight is favorable for larger birds and flapping flight for smaller birds, though the partitioning of time for each flight mode during migration is dependent on the tradeoff between time and energy (Hedenström, 1993; Duerr et al., 2015; Katzner et al., 2015; Miller et al., 2016). In theory, a timeminimizing migrator would be expected to fly with greater directional persistence and stronger directional bias than would an energy-expenditure minimizer. Such net energy maximizers would be expected to take advantage of en route foraging opportunities and may divert or delay to replenish energy reserves. (Note that "energy minimization" has been used to describe this strategy e.g., Alerstam, 2011; Miller et al., 2016, but we use "net energy maximization" for clarity.) If time is less important, a net energy maximizer is less restricted and can spend additional time seeking an energetically superior path; the emergent path would then be more tortuous with less directional bias toward the final destination at any given point along the route. Time minimization and net energy maximization strategies are not mutually exclusive, however, and the emergent strategy and behaviors in any given migrating individual lies along a continuum (Alerstam, 2011; Miller et al., 2016).

Obligate soaring migrants must also consider routes based on their energy landscape (Shamoun-Baranes et al., 2010), the energetic constraints of movement over space (Shepard et al., 2013), which also contributes to a migrant's location along the behavioral continuum. While soaring migrants can stopover, their energy landscape is more complex. Meteorological conditions are at least as important as foraging resources for soaring migrants, which can be extremely dynamic and subsidize the energetic cost of flight directly via uplift (Pennycuick, 1971; Alerstam, 1979; Spaar and Bruderer, 1997; Gill, 2007; Duerr et al., 2012; Murgatroyd et al., 2018).

The flight performance of soaring migrants relative to subsidies provided by meteorological conditions has been well documented (Pennycuick, 1971; Alerstam, 1979; Spaar and Bruderer, 1997), establishing a clear link between diurnal migrant behavior and development of the atmospheric boundary layer. Two primary forms of uplift arise by (1) wind interacting with topography to form upslope wind or mountain waves (air currents forming standing waves established on the lee side of mountains; hereafter orographic uplift) and (2) solar heating of the earth's surface to generate thermal uplift. Other forms arise from turbulent eddies over small landscape features and ocean waves modifying the air. The dynamic nature of atmosphericallydriven flight subsidies requires detailed movement data as well as carefully designed analytical techniques to investigate certain mechanisms hidden in those data.

Our understanding of migratory processes has advanced enormously in the past 30 years, as animal tracking technology developed from a novelty of coarse observation to a core method for observing animal behavior and movement in incredible detail (Luschi et al., 1998; Sawyer et al., 2005; Bridge et al., 2011; Katzner et al., 2015; Hooten et al., 2017). Global Positioning System (GPS) telemetry, in particular, allows remote observation of animal relocations across a broad spatiotemporal scale. GPS transmitters are now light and reliable enough to study the complete migrations of many large soaring migrants, including golden eagles Aquila chrysaetos, which often rely on flight subsidies during migration (Katzner et al., 2015). Golden eagles and other large soaring birds have been used as model systems for phenomenologically evaluating questions about migratory flight performance and migration strategies (sensu Duerr et al., 2012; Lanzone et al., 2012; Katzner et al., 2015; Vansteelant et al., 2015; Miller et al., 2016; Shamoun-Baranes et al., 2016; Rus et al., 2017). For example, Lanzone et al. (2012) and Katzner et al. (2015) found that golden eagles use both thermal and orographic uplift to subsidize migratory flight, although thermal soaring was often more efficient in long distance, directed flight (Duerr et al., 2012). While these studies have contributed to our understanding of soaring migration and have laid a foundation for more detailed approaches, they relate meteorology to derived movement metrics, rather than incorporate them into process-based models that mechanistically predict movement, and ignore the temporal dependence between serially observed locations (i.e., autocorrelation). Not accounting for such autocorrelation imparts bias on certain estimated parameters (e.g., variances) thereby affecting inference through, for example, underestimating uncertainty. Consequently, the links between resources distributed over the landscape, such as flight subsidies, and behavioral budgets, including stopover behavior, during migrations of soaring birds remain unclear.

Unlike previous approaches, process-based, mechanistic movement models allow explicit inference of the underlying mechanisms driving movement (e.g., changes in behavior) that may not be available from conventional phenomenological analytical approaches (Turchin, 1998; Nathan et al., 2008; Hooten et al., 2017). While it is impossible to understand fully the intricacies in animal movement, we can pose mathematical models (e.g., correlated random walks) to approximate the movement process (Kareiva and Shigesada, 1983; Turchin, 1998). We can then fit these models statistically to observed data to estimate parameters describing behavior and its relationship with dynamic environmental features that moving animals experience (Blackwell, 1997, 2003; Morales et al., 2004; Breed et al., 2017; Hooten et al., 2017). Many of the recently developed mechanistic movement models are built in a discrete stateswitching framework, where animals switch between discrete behavioral states (see Hooten et al., 2017, and references cited therein). Choosing both the biologically relevant and quantitatively supported number of states, as well as interpreting the identified states in a biological context, remains challenging (Patterson et al., 2017; Pohle et al., 2017). Often, this challenge leads researchers to artificially limit the number of states and/or collapse two or more states into one biologically interpretable state. For example Pirotta et al. (2018), presented a model with five discrete kinds of avian flight, but the complexity of the model made interpreting those states difficult and poorly matched classifications manually identified by an expert.

In many cases, a more natural approach to modeling an animal's movement process is along a dynamic continuum, rather than as switching between discrete behavioral states (Breed et al., 2012; Auger-Méthé et al., 2017; Jonsen et al., 2019). Modeling along a continuum may be an especially useful approach for understanding movement behavior in soaring birds, considering the dynamic nature of atmospheric processes that influence movements. Here, we developed and applied a flexible mechanistic movement model based on a correlated random walk with time-varying parameters. This novel model was fit to movement data collected via GPS telemetry to understand how individuals in a population of long-distance soaring migrants use flight subsidies and budget stopover and migration behavior. Specifically, we were interested in identifying which flight subsidies influence stopover and migratory behavior and how the effect of key subsides and behaviors varied between spring and fall migrations. Our approach resembled continuous-time correlated random walks (Johnson et al., 2008; Blackwell et al., 2015; Gurarie et al., 2017; Michelot and Blackwell, 2019), but was easily implemented and yielded a relatively small number of dynamic parameters that could be directly interpreted biologically. A set of candidate models could be ranked, with model selection approaches, providing inference on how behavioral budgets and meteorological variables interacted to give rise to the observed migration paths. Modeling the effects of dynamic wind and uplift variables as time-varying movement behaviors of migratory golden eagles further allowed new details to emerge without imposing artificially discrete states.

#### METHODS

#### Model System

The golden eagle is a large, soaring raptor, distributed across the Holarctic (Watson, 2010). Golden eagles are predatory and opportunistic, utilizing many taxa for food resources, ranging from small mammals and birds to ungulates, often scavenging carrion (Kochert et al., 2002; Watson, 2010). While many populations are classified as partial migrants, most individuals that summer and breed above approximately 55◦N in North America are considered true long-distance migrants (Kochert et al., 2002; Watson, 2010). The population we observed in this study migrates over the mountainous regions of western North America between a breeding range primarily in southcentral Alaska, USA and a broad overwintering range in western North America that ranges from the southwestern US to central British Columbia and Alberta, Canada (Bedrosian et al., 2018).

#### Data Collection

We captured golden eagles with a remote-fired net launcher, placed over carrion bait near Gunsight Mountain, Alaska (61.67◦N 147.35◦W). Captures occurred between mid-March and mid-April 2014-2016. Fifty-three adult and sub-adult eagles were equipped with 45-g back pack solar-powered Argos/GPS platform transmitter terminals (PTTs; Microwave Telemetry, Inc., Columbia, MD, USA). Eagles were sexed molecularly and aged by plumage.

PTTs were programmed to record GPS locations on duty cycles, ranging from 8 to 14 fixes per day during migration, depending on year of deployment. PTTs deployed in 2014 were set to record 13 locations at 1-h intervals centered around solar noon plus a location at midnight local time. PTTs deployed in 2015 were programmed to record 8 locations with 1-h intervals centered around solar noon, and PTTs deployed in 2016 took eight fixes daily at regular 3-h time intervals. Note that the PTTs deployed in 2015 did not record locations overnight. Poor battery voltage from September to March often resulted in PTTs failing to take all programmed fixes, so the resulting GPS tracks had missing observations during these periods. Tags lasted multiple seasons, and in fact many are still deployed and transmitting at this writing. We chose to limit this analysis to the migrations that occurred in 2016. The spring and fall migratory pathways of the 2016 migration from 26 tags were available and suitable for analysis in that year: 11 deployed in 2014, 7 deployed 2015, and 8 deployed 2016. Tracks were suitable for analysis based on having few missing data, with no more than a few days of consecutive missing locations.

Movement data were managed in the online repository Movebank (https://www.movebank.org/), and we used the Track Annotation Service (Dodge et al., 2013) to extract flight subsidy (wind and uplift) data, specific to each PTT location and time of recording that location, along eagle tracks. The Track Annotation Service derives uplift variables from elevation models and weather and atmospheric reanalyses (Bohrer et al., 2012). We followed the Movebank recommendations for interpolation methods; details are below.

#### Movement Model

We developed a correlated random walk (CRW) movement model to reveal how changes in behavior give rise to the movement paths of migrating eagles. We chose to use a dynamic, time-varying correlation parameter, which represents behavior as a continuum rather than discrete categories, to capture complex behavioral patterns that could occur on multiple temporal and spatial scales (Breed et al., 2012; Auger-Méthé et al., 2017; Jonsen et al., 2019). We believe this approach can offer substantial flexibility, as a continuous range of behaviors is more realistic and, as we show, more naturally allows modeling behavior as a function of covariates.

The basic form of the model was a first-difference CRW presented by Auger-Méthé et al. (2017), which can take the form:

$$\mathbf{x}\_{i}|\mathbf{x}\_{i-1}, \mathbf{x}\_{i-2} \sim \mathcal{N}\_{2} \left( \mathbf{x}\_{i-1} + \boldsymbol{\nu}\_{i} \frac{\Delta t\_{i}}{\Delta t\_{i-1}} (\mathbf{x}\_{i-1} - \mathbf{x}\_{i-2}), \ \boldsymbol{\Sigma}\_{i} \right), \tag{1}$$

where

$$
\Sigma\_i = \begin{bmatrix}
\Delta t\_i^2 \sigma\_\chi^2 & 0 \\
0 & \Delta t\_i^2 \sigma\_\chi^2
\end{bmatrix}, \quad \sigma\_\chi, \sigma\_\mathcal{\mathcal{V}} > 0. \tag{2}
$$

Here, 1ti = ti − ti−1 represents the time interval between Cartesian coordinate vectors **x**i and **x**i−1 for the observed locations of the animal at times ti and ti−1. Incorporating autocorrelation in behavior, γi constitutes a random walk, such that

$$
\gamma\_i | \gamma\_{i-1} \sim \mathcal{N} \left( \mathcal{Y}\_{i-1}, \ \Delta t\_i^2 \sigma\_\upsilon^2 \right), \ \ \sigma\_\upsilon \gg 0. \tag{3}
$$

γi correlates displacements (or "steps") and can be interpreted to understand the type of movement, and thus behavior, of migrating individuals: estimates of γi closer to one indicate directionally-persistent, larger-scale migratory movement, while estimates of γi closer to zero indicate more-tortuous, smallerscale stopover movement (Breed et al., 2012; Auger-Méthé et al., 2017). Scaling γi by 1ti 1ti−1 and the variance components by 1t 2 i allows us to accommodate unequal time intervals (Auger-Méthé et al., 2017), which can arise from a PTT's pre-programmed duty cycles and/or missed location attempts. This assumes that over longer time intervals an animal is likely to move greater distances and that the previous step will have less influence on the current step. Notably, in introducing 1ti , this CRW essentially becomes a correlated velocity model presented in terms of displacement vectors (**x**i−1−**x**i−2) (Johnson et al., 2008; Blackwell et al., 2015; Gurarie et al., 2017), most closely resembling the autocorrelated velocity model presented by Gurarie et al. (2017). Because location error of GPS data is negligible compared to the movement of most large vertebrates (Hooten et al., 2017), we did not incorporate an observation equation to handle location error. While a covariance parameter could be added to the model, we chose to fix covariance to zero (equation 2), which assumes that movement in the x and y dimensions are independent. This assumption has been suggested to be potentially problematic (Dunn and Gipson, 1977; Blackwell, 1997); however, it is common and has been shown to draw reasonable inference from real data, as well as recover known parameters from simulated data (Breed et al., 2012, 2017; Auger-Méthé et al., 2017; Jonsen et al., 2019). To support this, we compared results from the model assuming zero covariance to one fit assuming equal variance in x and y—like estimating covariance, this ensures invariance under linear transformation of the coordinate system—to illustrate that inference remains unaffected by this assumption (**Appendix 1**).

Extending this CRW to introduce environmental covariates, we first made the assumption that an individual's behavior can be adequately explained by the previous behavior plus some effects of environmental conditions and random noise. This modeling approach and philosophy aligns with the movement ecology paradigm presented by Nathan et al. (2008): An animal's movement path is influenced by its internal state and the environmental conditions it experiences. We modified the behavioral (or internal state) process—previously described above as a pure random walk in one dimension (Equation 3)—similar to a linear model with a logit link function. The logit link constrains γi ∈ [0, 1] and allowed us to model it as a linear combination of continuously-distributed random variables (Jonsen et al., 2019). These variables were different meteorological conditions affecting flight subsidies. Now,

$$\gamma\_i' = \log\left(\frac{\gamma\_i}{1-\gamma\_i}\right),\tag{4}$$

where

$$\boldsymbol{\gamma}\_{i}^{\prime} = \boldsymbol{\gamma}\_{i-1}^{\prime} + \mathbf{Z}\_{i}^{\mathrm{T}} \boldsymbol{\mathfrak{B}} + \boldsymbol{\epsilon}\_{i},\tag{5}$$

$$
\epsilon\_i \sim \mathcal{N}\left(\mathbf{0}, \,\,\Delta t\_i^2 \sigma\_\upsilon^2\right),
\tag{6}
$$

and **Z** T i is the row vector of environmental covariates associated with **x**i . Each element of the vector β is an estimated parameter representing the magnitude and direction of the effect of its respective covariate on the correlation parameter γi in addition to the effect of γi−1. Note that including γi−1 here preserves explicit serial correlation in the behavioral process so that any additional environmental effect is not overestimated. γ ′ i is only used to estimate γi ; any behavioral interpretations are made in terms of γi .

#### Model Fitting Subsetting Tracks

Of the 26 eagles producing suitable data in 2016, we fit the model to 15 spring and 16 fall adult golden eagle migration tracks recorded by 18 adult males and 8 adult females in 2016. This included both spring and fall migrations for five individuals. In reporting the results, we assumed any individual random effects of including both migrations for these few individuals to be negligible, which seems reasonable given fitted parameters presented in **Table S1** in Appendix. The model was fit only to the migratory periods, plus two fixes prior to departure to ensure valid parameter estimates at the onset of migration. Data were constrained to migratory periods under the following rules: The first migration step was identified as the first directed movement away from what was judged to be an individual's summer (or winter) range with no subsequent return to that range, and the final migration step was defined as the step terminating in the apparent winter (or summer) range. This assignment was usually straightforward; however, in some cases there were apparent pre-migration staging areas. These were not considered part of migration and excluded from the analysis here; movement data from these individuals collected during the breeding and overwintering periods are neither presented nor analyzed here.

#### Environmental Covariates

Golden eagles can switch between using thermal and orographic uplift as flight subsidies (Lanzone et al., 2012; Katzner et al., 2015), so we included both variables as covariates affecting the correlation parameter in the behavioral process of the CRW (equation 5). Thermal uplift **z**tu and orographic uplift **z**ou are measured in m/s with **z**tu, **z**ou ∈ [0,∞). Thermal uplift was bilinearly interpolated from European Center for Medium-Range Weather Forecasts (ECMWF) reanalyses, and orographic uplift from the nearest neighbor (grid cell) by pairing National Center for Environmental Predictions (NCEP) North American Regional Reanalysis (NARR) data with the Advanced Spaceborne Thermal Emission Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM; Brandes and Ombalski, 2004; Bohrer et al., 2012). We also introduced wind as a covariate in the behavioral process, as it can influence eagle flight as well as the flight and energy landscape of many birds during migration (Shamoun-Baranes et al., 2017). Wind data were bilinearly interpolated from the NCEP NARR u (easterly/zonal) and v (northerly/meridional) components of wind predicted 30 m above ground in m/s, from which we calculated the wind support **z**tw, such that **z**tw ∈ (−∞,∞) (Safi et al., 2013), where positive values correspond to tailwind and negative values headwind. The bearings used to calculate each ztw,i were the compass bearings required to arrive at **x**i+1 from **x**i .

We included a time of day interaction in the model because of clear diurnal effects. This also helped reduce zero inflation, particularly for thermal uplift, which often decays to zero after sunset due to heat flux and atmospheric boundary layer dynamics. To introduce the interaction, we used a dummy variable **z**0, such that z0,i = 0 when ti fell after sunset but before sunrise and z0,i = 1 when ti fell after sunrise but before sunset. This assumed behavior was not dependent on the covariates at night—the combination of covariates becomes zero when z0,i = 0—which is sensible given observed diurnal behavioral cycles. Sunrise and sunset times local to each GPS point were calculated in R with the "sunriset" function in the package "maptools" (Bivand and Lewin-Koh, 2016; R Core Team, 2016). Writing out the matrix operation in Equation (5), the final overall formulation of the behavioral process for the full model was:

$$\begin{split} \boldsymbol{\chi}\_{i}^{\prime} &= \boldsymbol{\chi}\_{i-1}^{\prime} + \left[ \beta\_{0} + \beta\_{\boldsymbol{\alpha}\boldsymbol{\epsilon}} (\boldsymbol{z}\_{\boldsymbol{\alpha}\boldsymbol{u},i} \times \boldsymbol{z}\_{0,i}) + \beta\_{\boldsymbol{\text{tu}}} (\boldsymbol{z}\_{\boldsymbol{\iota}\boldsymbol{u},i} \times \boldsymbol{z}\_{0,i}) \right. \\ &\left. + \ \beta\_{\boldsymbol{\text{tw}}} (\boldsymbol{z}\_{\boldsymbol{\iota}\boldsymbol{w},i} \times \boldsymbol{z}\_{0,i}) \right] + \boldsymbol{\epsilon}\_{i}, \end{split} \tag{7}$$

Prior to fitting the model, we followed Gelman et al. (2008) and log-transformed the uplift covariates and standardized variance to 0.25. We used a shifted log-transformation (Fox and Weisberg, 2019); adding one to the covariates prior to the logtransformation preserved zeros (i.e., zeros mapped to zero under the transformation). The distribution of raw wind support data appeared Gaussian, so it was only centered and standardized.

#### Parameter Estimation and Model Selection

We fit our correlated random walk (CRW) in a Bayesian framework. Because the model has explicit serially correlated parameters, we used Hamiltonian Monte Carlo (HMC) over more conventional Markov-chain Monte Carlo (MCMC; e.g., Metropolis steps) to sample efficiently from a posterior with such correlation.

Gelman et al. (2008) suggested Cauchy priors for logistic regression parameters; however, Ghosh et al. (2015) found that sampling from the posterior can be inefficient due to the fat tails of the Cauchy distribution. We thus chose Studentt priors centered on zero (µ = 0 and σ = 2.5) with five degrees of freedom as weakly informative priors for the covariate parameters. Weakly informative normal priors were placed on the variance parameters of the model.

We implemented HMC with R and Stan through the package "rstan" (R Core Team, 2016; Stan Development Team, 2016). Working R and Stan code, including details on prior choice, and example data are provided as **Supplementary Material**, as well as supplementary tables and figures (**Appendices 1–3**). The model was fit to each track independently with five chains of 300,000 HMC iterations, including a 200,000 iteration warm-up phase, and retaining every tenth sample. Convergence to the posterior distribution was checked with trace plots, effective sample sizes, posterior plots of parameters, and Gelman diagnostics (Rˆ) for each model fit.

We compared candidate models with leave-one-out crossvalidation approximated by Pareto-smoothed importance sampling (PSIS-LOO) in R with the package "loo" (Vehtari et al., 2016, 2017). The candidate models included possible combinations of environmental covariates plus a null CRW model without covariates. To limit model complexity and because we were interested in competing hypotheses about key predictors of behavior, we chose not to include interactions beyond time-of-day. We ranked the models by the expected log pointwise predictive density (elpd; i.e., out-of-sample predictive accuracy) transformed onto the deviance scale (looic; Vehtari et al., 2017), which created a measure on the same scale as common information criterion (e.g., AIC) and allowed applying the rules of more traditional information-theoretic model selection (e.g., Burnham and Anderson, 2004). The model with the lowest looic was considered the best fit to the data, but if other models were within two looic of the top model, each, including the top model, were considered equally supported by the data.

To understand how the predictive ability of the full model varied among tracks, we also computed a pointwise average of the elpd for each track (Vehtari et al., 2017). Normalizing by the sample size allowed comparing the out of sample predictive ability of the full model across individual migration tracks (**Table S1** in Appendix). The elpd (and looic), being sums, are otherwise dependent on the sample size for each model fit.

### RESULTS

#### Model Performance and Diagnostics

We fit eight candidate formulations of our CRW model to 31 migration tracks, equating to 248 total model fits. Chain mixing, Gelman diagnostics (Rˆ) close to one, and large effective sample sizes for all parameters indicated convergence to the posterior for most model fits. Posteriors of parameters appeared symmetric, also indicating the model was well behaved (**Figure S1** in Appendix). Across all migration tracks, the full model showed strong evidence of convergence, but for five tracks, we did not consider the null model converged to the posterior (e.g., Rˆ > 1.01). The five migrations for which the null model did not converge were not included from formal model selection.

#### Behavior During Migration

Median (interquartile range) departure and arrival dates were 5 March (4.5 d) and 27 March (6.4 d) in the spring and 29 September (11.7 d) and 16 November (15.5 d) in the fall. On TABLE 1 | Summary statistics of flight subsidies encountered by migrating golden eagles that summer in Alaska.

Variables were interpolated in space and time from weather reanalyses to eagle locations recorded by GPS telemetry. Units for all variables are m/s.

agrand mean across discrete GPS locations with individual migration tracks pooled.

average, eagles encountered similar orographic uplift in spring and fall but more intense thermal uplift and tailwind in the spring (**Table 1**).

The model revealed that eagles changed their behavior on multiple scales. First, there were very strong daily rhythms in behavior during migration, with birds migrating or moving more slowly and tortuously during the day and stopping at night (**Figures 1**, **2**). Explicitly including a time-of-day interaction could cause a daily rhythm to emerge as an artifact of model specification. However, accounting for serial correlation in behavior (Equation 5) limited that possibility. Additionally, prolonged periods of movement without an apparent daily rhythm suggest that, where daily rhythms are observed they are not a product of model specification (**Figure 2**). Second, there was some evidence of stopover-like behavior, but with individuals continually moving along the migration route while exhibiting less directional persistence in movement (**Figure 3**). The continuation along the migration route while in a stopoverlike state is highlighted by track segments extended over space associated with low and intermediate estimates of γi (blue/purple in **Figure 3**).

There was also a clear effect of season on movement patterns and behavior. Spring was characterized by straighter, more direct trajectories and punctuated by slower, more tortuous, stopoverlike movement; whereas, fall movements were much more tortuous overall and regular patterns in changes in movement rate and/or tortuosity less clear (**Figures 1**–**3**). The distributions of estimated γ values also clearly indicate that daytime movements were most frequently directed migratory moves in the spring; whereas, in the fall, the bimodal distribution indicates more equivalent partitioning between directed migratory moves and slower stopover type movement, with significant time spent exhibiting behaviors associated with intermediate tortuosity and movement rate (**Figure 3**).

#### Environmental Covariates

While there were differences in some environmental covariates between spring and fall (**Table 1**), parameter estimates from the full model (all covariates) indicate that there was little to no difference in effect of flight subsidies (i.e., wind and uplift) on behavior between spring and fall (**Figure 4**, **Table S1** in Appendix). Including environmental covariates in the behavioral process, though, improved model fit for almost all fitted migrations (**Table 2**). Positive coefficients on the thermal uplift covariate indicate that increasing thermal uplift resulted in more highly-correlated displacements, or migratory movements. Despite that, there were some migration bouts not associated with great thermal uplift (**Figures 1**, **2**). Coefficients close to zero for orographic uplift and wind support indicate that, in general, they were not strong drivers of directionallypersistent movements.

Based on the model selection, the best-fitting formulation of the environmental drivers of the behavioral process was variable across individuals. However, in almost all cases, some form of flight subsidy was used and there was little difference between the spring and fall seasons in the pattern of subsidy use (**Table 2**). The high variability across individuals (**Table S1** in Appendix) was likely due to differing weather patterns and thus subsidy sources encountered and/or used by each eagle as migrations were not synchronous (in time or space) across individuals. In addition, inter-individual variation was much larger than any difference attributable to demographic variables; we found no evidence that difference in sex or age explained patterns of flight subsidy use during migration. Note, though, that all eagles included in this analysis were in adult plumage, so strong age effects would not necessarily be expected.

Comparing the pointwise elpd across tracks revealed that the out of sample predictive ability of the full model varied among individuals (**Table S1** in Appendix). It also showed that predictive ability was greater for fall migrations than spring.

### DISCUSSION

Here, we develop and demonstrate how dynamic parameter CRW models fit to GPS data reveal the effects of variable flight subsides available along migration routes. Use of these subsidies gives rise to diverse patterns in the movement of a longdistance soaring migrant. Behavioral changes occur continuously as available subsidies shift over time and space. These key driving mechanisms underlie emergent movement paths, yet such processes are often hidden in the discrete satellite observations available. Our mechanistic modeling approach allowed linking of dynamic meteorology to changes in behavior, and those changes in behavior to the observed movement paths, revealing time series of behaviors more complex than individuals simply apportioning time between migration and stopover.

#### Model Performance

Incorporating time-varying parameters into movement models has been a relatively infrequently utilized approach (Breed et al., 2012; Auger-Méthé et al., 2017; Jonsen et al., 2019). Here we provide a case study for its utility and developed the approach for achieving practical biological inference about movement processes. Modeling the serial correlation in movement as a function of environmental covariates (equation 4), allowed simultaneous inference of behavior and the effect of environmental covariates on behavior from animal trajectories with regular and irregular duty cycles and containing missing observations. While other methods exist to handle missing data, the behavioral patterns we found would be more difficult to reveal with a state-switching movement model (e.g., hidden Markov models (HMMs); Michelot et al., 2016) because each

FIGURE 1 | Time series of behavior parameter γ from correlated random walk model with full behavioral process (orographic uplift, thermal uplift, and wind support as predictors) for two golden eagles during spring migration with PTTs reporting on different duty cycles. Upper panel is 13 hourly centered on solar noon plus one at midnight, and the lower panel is 8 hourly centered on solar noon. γ close to one reflect movements associated with migratory behavior, and γ close to zero stopover behavior. Points are times of observations, and lines are linear interpolations between points. Hue indicates intensity of thermal uplift, with yellow indicating greater and blue lower. Note the daily rhythm in behavior associated with intense thermal uplift, stopover periods of one or more days, and the intermediate periods suggesting fly-and-forage.

FIGURE 2 | Time series of behavior parameter γ from correlated random walk model with full behavioral process (orographic uplift, thermal uplift, and wind support as predictors) for three golden eagles during fall migration with PTTs reporting on different duty cycles. Upper panel is 13 hourly centered on solar noon plus one at midnight, middle panel is 8 hourly centered on solar noon, and lower panel is fixed 3-h interval. γ close to one reflect movements associated with migratory behavior, and γ close to zero stopover behavior. Points are times of observations, and lines are linear interpolations between points. Hue indicates intensity of thermal uplift, with yellow indicating greater thermal uplift and blue lower. Note the daily rhythm in behavior and extended stopovers as well as periods intermediate values suggesting fly-and-forage.

step would be forced into a discrete behavioral state from a set of usually 2–3 discrete states. Moreover, although hidden-state models have been introduced that have more than five discrete states (e.g., McClintock et al., 2012), these states can require ancillary data streams (e.g., accelerometry) to discriminate and remain extremely difficult to employ and interpret in practice

FIGURE 3 | Golden eagle migration trajectories (N = 15 spring and N = 16 fall). Hue indicates value of behavioral parameter γ estimated with the correlated random walk model with full behavioral process, including orographic uplift, thermal uplift, and wind support as predictors. Insets show the relative frequencies of estimates of γ assigned to the displacements between observed daytime GPS locations. γ close to one reflect movements associated with migratory behavior, and γ close to zero stopover behavior. Daily rhythms, revealed in Figures 1, 2, are not apparent here because the birds moved so little at night.

FIGURE 4 | Point estimates of environmental covariate effect parameters (βou, βtu, βtw) on golden eagle behavior and movements during migration (N = 15 spring and N = 16 fall). Estimates are from the correlated random walk model with full behavioral process, including orographic uplift, thermal uplift, and wind support as predictors.

(Patterson et al., 2017; Pohle et al., 2017). Finally, as HMMs include greater numbers of potential states, they tend to fit better than models with fewer states as judged by classical model selection approaches, such as AIC, even when additional states TABLE 2 | Number of golden eagle migration tracks recorded by GPS transmitters that each candidate formulation of the behavioral process in the correlated random walk model fit the best, according to approximate leave-one-out cross-validation (Table S1).

"therm" corresponds to thermal uplift, "oro" to orographic uplift, and "twind" to wind support.

a tally given to model with lowest information criterion (looic; Vehtari et al., 2016); if one or more models were within two looic of the top model, each was given a tally. boro + therm + twind.

are neither biologically meaningful nor sensible (Pohle et al., 2017). Implementing models with dynamic parameters that can be interpreted along a behavioral continuum seems a more natural approach for many animal movement questions.

#### Effects of Tag Programming

While our CRW model revealed the same trends across duty cycles and was generally robust to the different duty cycles (**Figures 1**, **2**), the most detail in daily behavioral rhythms was revealed in tracks with a fixed 3-hr time interval (lower panel in **Figure 2**), as it provided data throughout the 24-h day at regular intervals. The other duty cycles were initially chosen to minimize the risk of battery depletion overnight. Although generally robust, duty cycles did affect model fitting. HMC permitted Bayesian inference rather efficiently for our model, considering elevated correlation in the posterior of parameters due to the model formulation. Model fits typically took no more than a few hours, though tracks with much more than several hundred locations sometimes took longer. Preliminary fitting of our model with Stan and Template Model Builder (TMB; following Auger-Méthé et al., 2017) suggested that Maximum Likelihood estimation (when fit with TMB) tends to fail frequently when tag programming results in uneven temporal coverage of each day (e.g., our 2015 duty cycle), while Bayesian inference still provided sensible parameter estimates in most cases. Although the model presented herein and the model presented by Auger-Méthé et al. (2017) can make up for irregular time intervals between observations, they do have limitations. Breed et al. (2011) offer an in-depth discussion of tag programming and its effects on model fitting and inference.

### Flight Subsidies as Drivers Migration of Behavior

Thermal uplift is a flight subsidy dependent on daily atmospheric boundary layer dynamics, and it was clearly an important driver of the daily rhythm in eagle movement (**Figure 4**, **Table S1** in Appendix). Intense thermal uplift was often associated with the peaks in daily migration bouts (**Figure 1**). The larger magnitude of the thermal uplift effect, relative to orographic uplift, was somewhat surprising, as many individuals in our sample followed the Rocky Mountains, a large potential source of orographic uplift. Golden eagles are known to use orographic uplift as a flight subsidy while migrating through the Appalachian Mountains in eastern North America (Katzner et al., 2015). Much of the Appalachians, however, is characterized by long, unbroken, linearly-oriented ridges. Wind blowing over these ridges produces long stretches of predictable orographic uplift (Rus et al., 2017). The Rocky Mountains, by contrast, are far more rugged and nonuniform, and conditions that might produce suitable upslope winds and mountain waves, as well as strong tailwinds, likely also generate violent turbulence and could impede efficient migratory flight. Soaring raptors have been shown to use small-scale turbulence to achieve subsidized flight (Allen et al., 1996; Mallon et al., 2016); however, unpredictable, nonstationary violent turbulence, which can occur in large, highelevation mountain ranges (Ralph et al., 1997), could produce unfavorable migratory conditions. The large effect of thermal uplift, thus, could indicate that the Rocky Mountains, a spine that spans almost the entire migration corridor for this population, as well as some areas further west (Bedrosian et al., 2018), serves as a network of thermal streets for migrating eagles (Pennycuick, 1998). More explicitly, intense sun on south facing slopes would be expected to generate linear series of thermals that birds could glide between during both spring and fall migration. It is important to keep in mind that the migrants could capitalize on fine-scale, localized features of certain flight subsidies, like orographic uplift and tailwind, that may not have been captured by the interpolated meteorological data used in our analyses. However, model selection for models including those variables did indicate they explained some variance in eagle movement, which we discuss further below.

Despite meteorological conditions along migration paths that differed between spring and fall and a stark difference between behavioral budgets, our results showed no clear difference in the use of flight subsidies between the spring and fall seasons (**Figure 4**). This finding contrasts with season-specific effects of flight subsidies on golden eagle migration shown phenomenologically in eastern North America, where thermal uplift was shown to be the key subsidy in migratory performance during spring, while wind with some additional support from thermal uplift is most important in the fall (Duerr et al., 2015; Rus et al., 2017). Although our results indicate that eagles use similar flight subsidizing strategies in both seasons, consistent with the differences from the eastern population, the actual behaviors performed during spring and fall migrations differed considerably. In spring, eagles used subsidies to drive a migration that allows timely arrival on the breeding grounds, consistent with a time minimization strategy. In the fall, flight was subsidized to minimize net energy use, which emerged as a much more diverse behavioral repertoire during a slower fall migration (**Figure 3**; Miller et al., 2016). The more rapid and direct flight punctuated by bouts of tortuous, stopover-like movement in the spring (**Figure 3**), suggest eagles pause, refuel, and/or perhaps wait for better migration conditions. This suggests eagles may employ, at least in part, a net energy maximization strategy (Hedenström, 1993; Miller et al., 2016), despite the need for timely arrival on the breeding grounds to avoid fitness costs (Both and Visser, 2001).

The behavioral time series of spring migrations showed some evidence of individuals responding less to thermal uplift as latitude increased (**Figure 1**). This likely corresponded to a general decay in thermal uplift as individuals migrated northward (**Supplementary Material**, **Figure S2** in Appendix). Reduced thermal uplift availability would be expected at higher latitudes due to the larger amounts of remaining spring snowpack and lower solar angles. Thus, golden eagles, and likely other soaring birds, migrating to high latitudes may need to budget behaviors carefully between time minimization and net energy maximization during spring migration to best take advantage of the reduced flight subsidy from thermal uplift and mitigate the greater energy demands of flight at higher latitudes.

While our results show that thermal uplift is the most important flight subsidy for the majority of migrating eagles sampled, the model selection indicated orographic uplift and wind support improved out of sample predictive accuracy and explained some variance in eagle movement. Additionally, variability in top models across individuals (**Table 2**, **Table S1** in Appendix) suggests among-individual variance in flightsubsidizing strategy. Although some of this variability can be attributed to real individual differences in behavioral strategy, it is at least as likely that individuals encountered different subsidies en route and used the subsidies they had available, as migrations across our sample were not synchronous. Given that orographic uplift and wind support parameter estimates were negative or close to zero for many individuals (**Figure 4**, **Table S1** in Appendix), those covariates likely predicted the periods of slower, more tortuous movements (i.e., γi closer to zero). Wind support occurred in top models for more individuals in the fall (**Table 2**, **Table S1** in Appendix), which is consistent with findings from others (McIntyre et al., 2008; Rus et al., 2017) and suggests it may be important during southbound migrations. Additionally, although there was variance among the types and combinations of subsidies used, the null model (without flight subsidies) was the best fitting model for very few tracks (**Table 2**, **Table S1** in Appendix), evidence that weather and flight subsidies are of importance to migrating golden eagles, and likely also to the migrations of similar soaring species. Lastly, we found that the full model had, on average, better predictive ability in fall than spring (**Table S1** in Appendix), suggesting that the weather variables explained more of the variance within movement paths in fall compared to spring; during spring migration, other internal state variables associated with greater time limitation that were not explicitly accounted for in the models were likely responsible for this seasonal difference.

#### Daily Rhythm and Migratory Pace

The full movement model revealed two clear, nested behavioral patterns in the long-distance migrations of golden eagles. First, there was a daily rhythm where inferred directed migratory movements (i.e., γi close to one) occurred most frequently around midday or early afternoon (**Figures 1**, **2**). Mechanistic models of animal movement have revealed diel behavioral rhythms in other taxa (Jonsen et al., 2006). The basic aspects of daily rhythms in vertebrate behavior have hormone controls (Cassone, 1990), but the benefits can include balancing migration progress and foraging bouts (Newton, 2008). Soaring migrants also benefit by synchronizing diel movement patterns with diel atmospheric cycles. That is, consistent with our results, diurnal soaring migrants express a general circadian behavioral rhythm, where flight performance and behavior is strongly tied to thermal development of the planetary boundary layer to take best advantage of atmospherically generated flight subsides (Kerlinger et al., 1985; Leshem and Yom-Tov, 1989; Spaar and Bruderer, 1996, 1997; Mateos-Rodríguez and Liechti, 2012).

The second behavioral pattern revealed was a general stopover pattern, whereby eagles changed behavior for one to several days while en route (**Figures 1**–**3**). These changes were consistent with searching movements (i.e., γ intermediate or close to zero), possibly representing foraging behavior. In terms of soaring raptors, however, very few reports of movement patterns and behavior during stopovers have been published. Stopover segments have been previously identified by speed or some other metric calculated from tracks, then excluded from subsequent analyses (e.g., Katzner et al., 2015; Vansteelant et al., 2015); occasionally, authors noted apparent enhanced tortuosity but explored it no further (e.g., Vansteelant et al., 2017). On occasions where stopover behavior was considered, classifications based on stay duration and travel distance or speed with hard, often arbitrarily chosen cutoffs between migrating and stopover segments were used (Chevallier et al., 2011; Katzner et al., 2012; Duerr et al., 2015; Miller et al., 2016). In contrast, our modeling framework aligned with the movement ecology paradigm (Nathan et al., 2008); it used the observed data— GPS locations, rather than a derived metric—and a theoretical movement process to infer behavior from movement patterns along tracks on a spectrum ranging from stopped to rapid, directionally-persistent movement.

Our analyses, however, showed that eagles still tended to continue along their migration route during periods of movement most resembling stopover, but with reduced movement rate and directional persistence (**Figures 1**–**3**). This pattern suggests a joint migration/opportunistic foraging behavior that resembles fly-and-forage migration (Strandberg and Alerstam, 2007; Åkesson et al., 2012; Klaassen et al., 2017), which is consistent with observations of en route hunting behavior of golden eagles by Dekker (1985). Such behavior could be used to maintain balance between time expenses and energy intake, as it allows simultaneous migration progress and foraging.

This pattern does not fall very well within the "stopover" paradigm (Gill, 2007; Newton, 2008), however, as true stops during the migrations we observed were rare, except for expected nightly stops. Rather, migrants seemed to change their pace—either by slowing down, moving more tortuously, or both—but still generally moved toward their migratory destination (**Figures 1**–**3**). Thus, instead of a discrete behavioral framework, whereby migrants switch between two migratory phases (migration and stopover) with very different movement and behavioral properties, we propose that, for certain taxa, a continuous alternative framework "migratory pacing" may be more appropriate and a natural way to interpret en route migratory behavior and movement dynamics. Such taxa would include some and perhaps many soaring migrants, as well some migrating species in other fluid environments such as fishes and marine mammals. Soaring birds, even when energy reserves are relatively depleted, likely can still make steady progress toward a migratory goal when flight subsidies are available. Flapping migrants, on the other hand, would not be able to achieve this as readily, due in part to the greater energy demands for sustained flight, and would require more regular refueling stopovers where migration progress is temporarily completely arrested. Both opportunism in foraging and use of energetic subsidies are likely key characters of fly-and-forage behavior and the ability to change pace of migration without actually stopping, as they relax the need for individuals to stopover in specific, food-rich habitats, which are required by most migrants with less flexibility in food and that lack the morphological specialization to maximally exploit the energetic subsidies available in moving fluids (Gill, 2007; Piersma, 2007).

Our model results revealed seasonal variability in migratory pacing by golden eagles. The tendency for eagles to exhibit movements matching fly-and-forage behavior, and pace their migrations more slowly was most apparent during fall migration. In contrast, spring migration was usually composed of much more punctuated events of slower-paced movements but these were still extended over space (**Figure 3**), indicating the eagles pace their migration and employ a mixed behavioral strategy to some extent in spring as well. During spring, hibernating mammalian prey would be minimally available, leaving carrion, along with a few non-migratory and -hibernating species (e.g., ptarmigan Lagopus spp. and hare Lepus spp.), as major food sources, which could help explain the more punctuated bouts of slower-pacing. Alternatively, individuals could have been slowed by poor weather conditions (Rus et al., 2017). Scavenging large ungulate carcasses would be extremely rewarding in terms of energy accumulation. Much of the carrion we used successfully to capture eagles was large ungulate (e.g., moose Alces alces), strongly suggesting that the population we sampled uses carcasses during migration. The bimodal distribution for the behavioral parameter γ in fall shows that eagles tended to budget daytime behaviors approximately equally between rapid, directed and slower-paced movements (**Figure 3**); the high frequency and range of intermediate values are, again, evidence for the more complex fly-and-forage and pacing dynamic, rather than eagles simply switching between stopover and migration. This behavioral complexity might be biologically important, allowing eagles to arrive on winter home ranges in better condition compared to migration strategies that do not incorporate en route foraging opportunity. In contrast to fall, daytime movements in the spring were typically faster-paced (i.e., larger-scale and directionally-persistent; **Figure 3**), consistent with a time minimization strategy, where eagles need to partition time more in favor of migration progress to ensure timely arrival on breeding grounds (Hedenström, 1993; Alerstam, 2011; Miller et al., 2016). We thus see in eagles, and propose more generally, that such pacing varies between and within seasons along the continuum between time minimization and net energy maximization strategies (Alerstam, 2011; Miller et al., 2016). A migrant's pace would be expected to depend upon their energetic demands, energetic subsidies available from the environment, and the importance of arriving at the migration terminus in a timely fashion (Nathan et al., 2008).

#### Implications and Conclusions

We developed and applied a movement model with timevarying parameters to help reveal the mechanisms underlying the migration of a long-distance soaring migrant that relies on incredibly dynamic flight subsidies. We found that variation in flight subsidies gives rise to changes in migrant behavior with thermal uplift seemingly most important. While these findings might be expected given previous phenomenological analyses (e.g., Duerr et al., 2012; Lanzone et al., 2012; Katzner et al., 2015; Vansteelant et al., 2015; Miller et al., 2016; Shamoun-Baranes et al., 2016; Rus et al., 2017), we were able to show how meteorology is a mechanism influencing changes in movement patterns and thus behavior.

In the behavioral budgets of migrating golden eagles, we identified an expected daily rhythm, as well as evidence for behavioral dynamics that would allow nearly simultaneous foraging and migration, which is greater complexity than the traditional stopover paradigm allows. Migratory pacing, facilitated by fly-and-forage behavior, expands the traditional notion of stopover, whereby a bird migrates until resting and refueling is required, at which point it stops for a brief period in specific habitat suitable for efficient foraging (Gill, 2007; Newton, 2008). This advance was enabled by incorporating time-varying parameters into the movement model, which revealed new behavioral patterns during migration of long-distance soaring migrants. While time-varying, dynamic parameters have been infrequently employed in movement modeling (Breed et al., 2012; Jonsen et al., 2013; Auger-Méthé et al., 2017), we have shown it is a promising approach that can overcome certain limitations in discrete state-switching models and help provide novel insight into animal behavior.

This approach also has potential for further development and for revealing additional new patterns in soaring bird movement; it has already been shown to help provide new insight for other taxa as well (Jonsen et al., 2019). Although we demonstrated the approach for several individual eagles, applying our methods across a larger sample and across more years will increase the inferential strength of our results. For example, previous work found effects of wind support and orographic uplift (e.g., Katzner et al., 2015; Vansteelant et al., 2015), where we, in accounting for an eagles' underlying movement process and the inherent autocorrelation in that process, found that those meteorological variables may be of less importance, at least compared to thermal uplift. It remains unclear though, whether these are systemspecific findings or a more general result. Additionally, the model we present has potential to help assess effects of habitat on the movement decisions of soaring birds and other species. One potential avenue for such would be incorporating the model into a resource selection framework (e.g., step selection function). Furthermore, given the movement process is parameterized in terms of coordinate vectors, the position likelihood could be straightforwardly extended to include the z axis to investigate questions regarding flight height of soaring birds or dive depth of marine species, assuming data of acceptable temporal resolution and location error are, or become, available.

### DATA AVAILABILITY

All movement data used for this manuscript are managed in the online repository Movebank (https://www.movebank.org/; IDs 17680093 and 19389828). The data contain information considered sensitive by the State of Alaska, but they could be made available at the discretion of the Alaska Department of Fish & Game and U.S. Fish & Wildlife Service. Code to fit the movement model to data and example data are provided as **Supplementary Material**.

#### ETHICS STATEMENT

Field procedures were conducted following the ADF&G Animal Care & Use Committee protocol #2013-036 and University of Alaska Fairbanks Institutional Animal Care & Use Committee protocol #859448.

### AUTHOR CONTRIBUTIONS

JE and GB conceived the ideas of the research presented herein. JE, MA-M, and GB designed analytical methods. TB, CB, and SL designed the field methodology. TB, CB, SL, and JE collected the data. JE analyzed the data, and led the manuscript. All authors contributed to drafts and editing of the manuscript and provided final approval for publication.

#### ACKNOWLEDGMENTS

We thank M. Kohan, B. Robinson, T. & D. Hawkins and many others for support in the field and J. Liguori and N. Paprocki for help aging eagles. We also thank T. Avgar, P. Doak, T. Katzner, K. Kielland, C. McIntyre, A. Scharf, and P. Smouse for helpful comments on drafts of the manuscript. Funding was provided by the Alaska Department of Fish & Game (ADF&G) through the federal State Wildlife Grant Program, and the U.S. Fish & Wildlife Service (USFWS) provided PTTs

#### REFERENCES

and data. JE was supported by the Calvin J. Lensink Fund during part of the project. MA-M acknowledges the support of the Natural Sciences and Engineering Research Council of Canada. The findings and conclusions of this paper are those of the authors and do not necessarily represent the views of the USFWS. This manuscript has been released as a pre-print at Eisaguirre et al. (2019).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00317/full#supplementary-material

behavioral dynamics from animal tracks. Ecol. Model. 235–236, 49–58. doi: 10.1016/j.ecolmodel.2012.03.021

unit of animal movement: synthesis and applications. Movem. Ecol. 5, 1–18. doi: 10.1186/s40462-017-0103-3

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Eisaguirre, Auger-Méthé, Barger, Lewis, Booms and Breed. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Scales of Blue and Fin Whale Feeding Behavior off California, USA, With Implications for Prey Patchiness

#### Ladd M. Irvine\*, Daniel M. Palacios, Barbara A. Lagerquist and Bruce R. Mate

*Department of Fisheries and Wildlife, Marine Mammal Institute, Oregon State University, Newport, OR, United States*

#### Edited by:

*Stan Boutin, University of Alberta, Canada*

#### Reviewed by:

*Inger Suzanne Prange, Appalachian Wildlife Research Institute (United States), United States Kimberly Goetz, National Institute of Water and Atmospheric Research (NIWA), New Zealand*

> \*Correspondence: *Ladd M. Irvine<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution*

> Received: *02 February 2019* Accepted: *23 August 2019* Published: *10 September 2019*

#### Citation:

*Irvine LM, Palacios DM, Lagerquist BA and Mate BR (2019) Scales of Blue and Fin Whale Feeding Behavior off California, USA, With Implications for Prey Patchiness. Front. Ecol. Evol. 7:338. doi: 10.3389/fevo.2019.00338* Intermediate-duration archival tags were attached to eight blue whales (*Balaenoptera musculus*; four females, three males, one of unknown sex), and five fin whales (*B. physalus*; two females, one male, two of unknown sex) off southern California, USA, in summer 2014 and 2015. Tags logged 1-Hz data from tri-axial accelerometers, magnetometers, and a depth sensor, while acquiring Fastloc GPS locations. Tag attachment duration ranged from 18.3 to 28.9 d for blue whales and 4.9–16.0 d for fin whales, recording 1,030–4,603 dives and 95–3,338 GPS locations per whale across both species. Feeding lunges (identified from accelerometer data) were used to characterize "feeding bouts" (i.e., sequences of feeding dives with <60 min of consecutive non-feeding dives), within-bout behavior, and to examine the spatial distribution of feeding effort. Whales fed near the tagging locations (Point Mugu and San Miguel Island) for up to 7 d before dispersing as far south as Ensenada, Mexico, and north to Cape Mendocino, California. Dispersal within southern California waters differed by sex in both species with males undertaking offshore, circuitous excursions, while females remained more coastal, suggesting that movement patterns on the feeding grounds may not be exclusively related to energy gain. Feeding bout characteristics were similar for both species, with the median bout having 24 dives and lasting 3.3 h for blue whales (*n* = 242), and 19 dives while lasting 2.7 h for fin whales (*n* = 59). Bout duration was positively correlated with the number of feeding lunges per dive within a bout for both species, suggesting whales left poor-quality prey patches quickly but fed intensively for up to 34.9 h when prey was abundant. Feeding bouts occurred further apart as the distance from shore increased, but there was no corresponding difference in the number of feeding lunges per dive, suggesting the whales were feeding at the same rate throughout their range, but that prey was more dispersed in offshore waters. This may be evidence of two feeding strategies, with spatially aggregated foraging around highly localized, topographically forced upwelling centers nearshore, and more dispersed foraging in larger areas of elevated, but patchy, productivity offshore.

Keywords: data loggers, accelerometers, GPS tracking, blue whale, fin whale, feeding behavior, sexual segregation

### INTRODUCTION

Blue (Balaenoptera musculus) and fin (B. physalus) whales are the two largest species of cetaceans and both occur off the west coast of the United States (USA). Blue whales arrive seasonally beginning in the late spring-early summer and feed on aggregations of krill (Thysanoessa spinifera and Euphausia pacifica; Fiedler et al., 1998; Croll et al., 2005; Nickels et al., 2018) until the late fall-early winter when they migrate south to breeding areas off Baja California, Mexico, in the Gulf of California, and near an offshore oceanographic feature called the Costa Rica Dome (Mate et al., 1999; Bailey et al., 2010; Irvine et al., 2014). Although the seasonal movements of fin whales are less well-understood, they occur year-round off southern California (Stafford et al., 2009; Sirovic et al., 2015; Scales et al., 2017) suggesting they do not follow the typical baleen whale pattern of migrating to lower latitudes during the winter to breed and calve (see also Edwards et al., 2015; Geijer et al., 2016; Jiménez López et al., 2019). Fin whales primarily feed on the same species of krill as blue whales, but can also feed on small fish and squid (Pauly et al., 1998). Both species are listed in the USA as "Endangered" under the Endangered Species Act and consequently are labeled as "Depleted" and "Strategic" stocks under the Marine Mammal Protection Act. Additionally, blue whales are considered "Endangered" and fin whales "Vulnerable" according to the IUCN Red List of Threatened Species (Cooke, 2018a,b).

As with other rorquals, blue and fin whales feed by engulfing large volumes of water and schooling prey, then expelling the water through fibrous baleen plates to retain the prey (Goldbogen et al., 2017). This engulfment, termed "lunge feeding," occurs during a rapid acceleration through/into a school of prey and can happen multiple times per dive. The evolution of this feeding behavior is closely tied with the animal's large size, as it allows them to efficiently exploit highly concentrated prey patches (Goldbogen et al., 2011, 2012; Goldbogen and Madsen, 2018). In blue whales this feeding behavior has been shown to vary based on local prey concentrations, indicating they are able to adapt foraging strategies to maximize prey capture (Goldbogen et al., 2015) and energetic efficiency (Hazen et al., 2015).

While sexual differences in distribution and behavior are welldocumented in toothed whales (Bigg et al., 1990; Whitehead, 2003; Parsons et al., 2009), such dynamics are less wellunderstood for rorqual whales, which are largely non-sexually dimorphic, but are known to occur (e.g., Laidre et al., 2009). Sex-based differences in acoustic behavior have been documented for both blue and fin whales, with most of the more complex vocal repertoire of each species being produced exclusively by males and believed to be related to reproduction, while singular calls are produced by both sexes in the context of foraging (Croll et al., 2002; Oleson et al., 2007; Stimpert et al., 2015; Lewis et al., 2018). Further, spatial variability in the type and rate of calling has been documented in blue whales off southern California, suggesting there may be spatial separation of behavioral activities (Lewis and Sirovic, 2018), although a lack of direct observation makes the underlying process unclear, as both feeding and reproductive calls are recorded throughout the feeding season in southern California (Lewis and Sirovic, 2018). Considering that animals in southern California waters are exposed to a variety of regional stressors like ship strikes (Berman-Kowalewski et al., 2010; Redfern et al., 2013) and anthropogenic sound (Goldbogen et al., 2013; DeRuiter et al., 2017), a potential spatial segregation arising from behavioral differences between sexes could lead to disproportionate impacts (Sprogis et al., 2016) over short temporal scales (Pirotta et al., 2018) or on long-term population fitness (Pirotta et al., 2019).

To understand broad-scale ecological patterns arising from observed species distributions, including interactions among sympatric species, it is often relevant to understand how finescale behavior influences the larger trend (Evans, 2012). The movements and distribution of blue and fin whales at the scale of the eastern North Pacific Ocean have been described using satellite telemetry (Bailey et al., 2010; Irvine et al., 2014; Scales et al., 2017), ship based surveys (Redfern et al., 2013), and acoustic data (Burtenshaw et al., 2004; Sirovic et al., 2015). However, there is little data to explain how the broad-scale distributions of these animals are influenced by variations in behavior at more local scales. Here we used intermediateduration data-logging tags to examine the fine-scale movements and diving behavior of blue and fin whales off southern and central California over periods of multiple weeks. Our primary goal was to characterize how feeding behavior varies across this region, including potential differences between the sexes and between the two species, which occur in sympatry in this region. Our results indicate a strong spatial and temporal heterogeneity in blue and fin whale foraging behavior, with implications for prey patchiness and quality. As such, this study provides new insights into the underlying drivers of broad-scale movement and occurrence during the feeding season for two of earth's largest predators.

#### MATERIALS AND METHODS

#### Tag Configuration and Deployment

This study used Advanced Dive Behavior (ADB) tags, a configuration of the Wildlife Computers (Seattle, WA, USA) MK-10 time-depth recorder platform, which can record depth, temperature, and tri-axial accelerometer, and magnetometer data at 1-Hz resolution for multiple weeks before releasing from the whale for recovery (Mate et al., 2017). A Fastloc R GPS receiver and patch antenna were included in each tag to acquire GPSquality locations (Bryant, 2007), along with an Argos Platform Terminal Transmitter for sending GPS location and summarized dive data messages via the satellite-based Argos Data Collection and Location System. Complete details of ADB tag construction and design configuration are provided in Mate et al. (2017).

Tags were deployed on blue and fin whales in southern California waters during August 2014 and July 2015. Field work was supported by the 25.6-m research vessel R/V Pacific Storm. Tags were mounted in a semi-implantable stainless steel housing and deployed at close range (2–4 m) from a 6.7-m rigidhulled inflatable boat using the Air Rocket Transmitter System, a modified compressed-air line-throwing gun (Heide-Jorgensen et al., 2001; Mate et al., 2007). Tags were deployed 1–4.5 m forward of the dorsal fin of the whale and no more than 20 cm down from the mid-line, following the protocol described in Mate et al. (2007). Tags were programmed to release from their housing for recovery 21 d after the start of field operations or if the tag recorded no change in depth for 24 h, indicating the housing had been shed from the whale and sunk to the bottom with the tag still attached. Following release, tags were located and recovered using an uplink receiver that was capable of acquiring, decoding, and solving Argos-transmitted Fastloc GPS location messages sent by the tags in real time, along with information on the tag's general location and the rate and direction of drift. A discussion of factors affecting attachment duration and tag recovery is presented in Mate et al. (2017).

Whenever possible, skin and blubber samples were collected simultaneously to tagging (on the same surfacing) using a crossbow. Crossbow bolts were fitted with circular cutting tips 4 cm in length and 8 mm in diameter, which removed a plug of skin and blubber from the whale. Sex determination was made from skin samples by amplification of regions on the X and Y chromosomes (Aasen and Medrano, 1990; Gilson et al., 1998).

#### Data Collection and Transmission

Tags were programmed to collect data from the various sensors (depth, accelerometer, magnetometer, and temperature) at 1 Hz for the duration of the deployment. While deployed, acquisition of a Fastloc GPS location was attempted every 7 min or the next time the whale surfaced after this time had elapsed. All data were stored in an onboard archive for download after tag recovery. The dive summary data from the Argos transmissions were not used in this study, although locations estimated from Argos transmissions (Argos, 2016) were used to examine movements of whales whose tags were not recovered. Further details about ADB tag dive summary messages are described in Mate et al. (2017).

#### Data Analysis

Maps were made using ArcGIS R software ArcMapTM 10.3 by Esri and the ArcGIS Online Ocean Basemap. Data manipulation and analysis for this study was conducted using Matlab (The Mathworks Inc, 2015) and R (R Core Team, 2018). The same suite of analyses were used for both blue and fin whales, except as noted for cases when data for one species were too limited for complete analyses. To detect feeding lunges, the change in the acceleration vector ("Jerk") was calculated from the accelerometer data as the norm of the difference between consecutive acceleration values (Simon et al., 2012; Allen et al., 2016). The Jerk is thus a measure of rapid changes in acceleration and orientation of the tagged whale. Lunge-feeding events in rorquals are characterized by a peak in Jerk with a coincident increase in the roll angle for multiple seconds, as the whale typically accelerates and rolls when opening its mouth to engulf prey (Goldbogen et al., 2006; Simon et al., 2012; Allen et al., 2016). A subsequent minimum in the Jerk value, as the whale ceases most movement to expel the water and filter out prey, signals the end of the lunge. Together, these three criteria (Jerk maximum, increase in roll, and subsequent Jerk minimum) were used to identify feeding lunges in the data records based on the lunge detection methodology described by Allen et al. (2016), but modified to use a more conservative Jerk peak threshold and without the use of acoustic data. The large number of dives recorded by tags during their deployment periods precluded direct confirmation of each feeding lunge identified by the lunge detection algorithm for the entire data record. Instead, we implemented a validation protocol by randomly selecting 10% of dives from each track and visually reviewing them to estimate the percentage of correctly detected lunges (true positives), falsely detected lunges (false positives), and correctly identified feeding dives (dives with at least one true positive). Validation statistics were summarized for each tag and are presented briefly in the Results section, with a more detailed description in the **Supplementary Material**.

Dive summaries were generated for each track by isolating any submergence >10 m in depth (hereafter a "dive") from the tag's depth record and calculating maximum dive depth, dive duration, and the number of lunges that occurred during the dive. The dive end times were then matched to the nearest GPS location recorded by the tag, as locations were generally collected as the whale surfaced from a dive. If there was not a location within 10 min of a dive, a location was estimated by linear interpolation between the temporally closest GPS locations before and after the dive using the dive time to determine where on the line the location should fall. For tracks with less frequent GPS locations this resulted in linear segments of interpolated dive locations that do not represent the exact movement of the whale.

In order to distinguish between series of related feeding dives, a log-survivorship analysis (Holford, 1980; Gentry and Kooyman, 1986) was conducted on the time between feeding dives (i.e., dives with at least one detected lunge) for each tag record. In this case, the log-survivorship analysis graphically showed the number of feeding dive sequences (on a log scale) as a function of the time between them, and the goal was to identify a point along the curve where the number of feeding dive sequences stabilized as time between feeding dives continued to increase, suggesting a natural break in the data. Sequences of dives defined by this criterion were isolated and labeled "feeding bouts." To assess the horizontal extent of each feeding bout and the overall spatial scale of foraging effort, minimum convex polygons were created using the corresponding dive locations for bouts with at least three GPS locations (i.e., not including bouts with interpolated dive locations; see **Supplementary Figure 1**). We report feeding bout summary statistics including bout duration, time, and distance between consecutive bouts, and distance to the closest bout from the entire track. We also report summaries of dives within each feeding bout, including number of dives per bout, mean maximum dive depth and duration, mean number of lunges per dive, and the number of dives without a lunge. The univariate distributions of feeding bout metrics were assessed graphically using probability density plots, while inter-species comparisons of the distributions were made using Bhattacharyya's similarity coefficient (Bhattacharyya, 1943; Guillerme and Cooper, 2016). When Bhattacharyya's similarity coefficient between two distributions is <0.05, the distributions are significantly different, and when the coefficient is >0.95, the distributions are significantly similar. Values between these two thresholds can be used to indicate the probability of overlap between the distributions but cannot be used to determine if they are significantly different or similar (Guillerme and Cooper, 2016).

Trends among these metrics were formally assessed using generalized linear mixed models (GLMMs), with the tag number as a random-effect grouping variable to account for differences between individuals, a fixed-effect "species" indicator variable, and an interaction term between the predictor of interest and the species indicator variable when inter-species comparisons were desired (Bolker et al., 2009; Harrison et al., 2018). Variables were log-transformed as needed based on graphical assessment of the data. Models were implemented in R using the lme4 package (Bates et al., 2015) or GLMMadaptive (Rizopoulos, 2019) using a two part/hurdle model if it was necessary to account for zero inflated data (GLMMzi). P-values were derived using the lmerTest package (Kuznetsova et al., 2017). We report and discuss p-values from these GLMMs in the context of levels of support for the outcome (rather than as a binary threshold of significance), with p-values <0.01 offering strong support, pvalues between 0.01 and 0.1 offering suggestive, but inconclusive support, and p-values >0.1 offering no support (Gerrodette, 2011; Greenland et al., 2016; Wasserstein and Lazar, 2016).

To further examine fine-scale differences in feeding behavior between individuals, we isolated sections of tracks where two whales were in proximity to each other. Whales were considered to be in proximity when locations from one track were no more than 1 km away from another whale within 30 min from the time of the location. We present a representative map of these track sections and also report dive summary statistics for each section of track in close proximity to another whale for comparison between individuals.

#### RESULTS

During this study, eight blue whales were tracked for a median of 22.4 d (range = 18.3–28.9 d; **Table 1**) and five fin whales were tracked for a median of 14.2 d (range = 4.9–16.0 d; **Table 1**). Tags were deployed near Point Mugu in 2014 and off the west end of San Miguel Island in 2015 with the exception of blue Whale # 2015\_838, which was deployed near Point Mugu in 2015 (**Figure 1**). All eight blue whale tags and three of five fin whale tags were recovered. In two cases the tags were found on the shore several years after having been thought to be lost.

Tagged animals of both species dispersed as far north as Cape Mendocino in northern California, and as far south as Ensenada, Baja California, Mexico, during their tracking periods (**Figures 2**, **3** and **Supplementary Figure 4**). However, most of the tracks, and the majority of feeding dives (89%) occurred within southern California waters. Two blue whales (Whale # 2014\_5650 and Whale # 2015\_4177) and one fin whale (Whale # 2014\_5685) made a clockwise loop across a large portion of southern California waters. All three of these whales were male, while females remained closer to shore until leaving southern California waters (**Figures 2**, **3**).

The eight blue whale ADB tags recorded a median of 126 dives > 10 m in depth per day (range = 73–207 dives/d; **Table 1**), whereas the three recovered fin whale tags recorded a median of 106 dives/d (range = 84–210 dives/d; **Table 1**). The number of Fastloc GPS locations recorded by the tags varied widely, with recovered blue whale tags recording a median of 72 locations per day (range = 10–139; **Table 1**) and recovered fin whale tags recording a median of 46.5 locations per day (range = 1–99; **Table 1**). Some variability in the number of locations was likely due to different hardware configurations within the tags, as all tags using newer Fastloc v. 3 technology recorded at least 54 locations per day while only one tag using Fastloc v. 1 recorded more than 47 locations per day (**Table 1**; see also Mate et al., 2017).

Validation of the feeding lunge detection algorithm indicated a mean true positive rate of 70.5% (sd = 20.1%) and a mean false positive rate of 13.6% (sd = 10.0%). The mean percentage of correctly identified feeding dives (dives with at least one lunge) was 85.3% (sd = 16.8%), suggesting mis-classified lunges (both false positive or false negative) often occurred during dives with other correctly identified lunges. Additional details of the validation methodology and results are presented in the **Supplementary Materials**.

Tagged whales of both species generally made deeper dives during the daytime than at night (**Supplementary Figures 2, 3** and **Supplementary Table 1**), although there was high variability within and between individuals, and daytime surface feeding was recorded on multiple occasions both visually while in the field and in the data record. Daytime dive depths were deepest near San Miguel Island where maximum dive depths reached 362 m for blue whales and 365 m for fin whales. Dive durations were as long as 30.7 min for blue whales and 23.1 min for fin whales (**Supplementary Figures 2, 3** and **Supplementary Table 1**). Almost no feeding lunges were recorded during dives >20 min in duration for either species. Feeding activity (as measured by lunge-feeding events) generally took place during daylight hours, although nighttime lunges were recorded on some occasions for both species (**Supplementary Figures 2, 3** and **Supplementary Table 1**). Most blue whale feeding effort was concentrated near the tagging areas, with additional areas of elevated feeding effort southwest of San Miguel Island and offshore of central California (**Figure 2** and **Supplementary Figure 4**). Whales tagged near Point Mugu also heavily used the nearshore waters extending south to San Diego. Tagged fin whales showed a similar trend (**Figure 3** and **Supplementary Figure 4**).

Log-survivorship curves for all tags stabilized at 60 min or less (**Supplementary Figure 5**), so a criterion of at least 60 min with no feeding dives was used to differentiate between feeding bouts. A total of 242 blue whale feeding bouts and 59 fin whale feeding bouts were identified in the tag records (**Tables 2**, **3**). For blue whales, the median number of feeding bouts made per whale was higher in 2014 (median = 35, range = 22–38, 4 tags) compared to 2015 (median = 23, range = 17–45, 4 tags), despite tags remaining attached for a median of over 7 d longer in 2015 (**Table 1**). The number of fin whale feeding bouts was more similar across years (range = 13–25 bouts, 3 tags), although the more limited number and duration of tracks was too small for inter-annual comparisons. Across both years, feeding bouts were separated by a median of 5.7 h (range = 1–231.9 h) for

TABLE 1 | Deployment summary for ADB tags attached to eight blue and five fin whales in southern California waters during August 2014 and July 2015.

\**Data were transmitted through Service Argos, Inc.*

*Some tags were not recovered due to poor weather and the tag's distance to shore after release (*>*160 km in some cases) limiting recovery opportunities. See Mate et al. (2017) for additional details.*

Median (recovered tags)

blue whales and a median of 3.1 h (range = 1–227.6 h) for fin whales, and were generally small in area (median = 5.6 km2 , range = 0.003–546.5 km2 for blue whales; and median = 7.6 km2 , range = 0.001–317.1 km2 for fin whales). The median blue whale feeding bout contained 23.5 feeding dives over 3.3 h (range = 4– 360 feeding dives and 0.2–34.9 h, respectively), while the median fin whale feeding bout contained 19 feeding dives and lasted 2.7 h (range = 4–142 feeding dives and 0.3–19.6 h, respectively; **Tables 2**, **3**).

Median bout duration for the four female blue whale tags (3.9 h) was twice that of males (1.8 h) although median bout duration of Whale # 2014\_5644 (a female) was = 1.9 h. Female feeding bouts had a lower proportion of non-feeding dives than those of males (0.26 vs. 0.38; **Table 2**). The distribution of feeding bout duration was similar for both blue and fin whales (Bhattacharyya's similarity coefficient = 0.934), with a strong peak near 2 h and a secondary peak at 14–15 h (**Figure 4**). Mean feeding lunges per dive within bouts varied substantially for both blue and fin whales, but feeding bout duration increased with increasing mean lunges per dive (GLMM, pvalue <0.001), with no significant difference between blue and fin whales (GLMM, interaction p-value = 0.84; **Figure 5**). The median bout duration for one blue whale (Whale # 2015\_840, of unknown sex) was 12.2 h, suggesting it fed almost continuously during daylight hours on many days. Another blue whale (# 2015\_838) fed continuously for almost 1.5 days (**Table 2** and **Figure 2**).

Comparison of the distributions for closest feeding bout and for the distance to shore of a feeding bout between blue and fin whales (**Figure 4**) using Bhattacharyya's similarity coefficient was inconclusive in terms of providing support for similarity or difference between them, although the probability of overlap was relatively high in both cases (Bhattacharyya's similarity coefficient = 0.840 and 0.833, respectively). Investigation of the distance between closest feeding bouts required accounting for zero-inflated data, as 44% of bouts overlapped spatially with at least one other bout from the same track, resulting in many zero distances. As the distance to shore of a feeding bout increased, the distance between closest feeding bouts also increased (GLMMzi, p-value <0.005), indicating that feeding bouts were more dispersed further offshore (**Figure 5**). However, there was no significant difference in the average number of lunges per dive made within bouts as a function of the distance from shore (GLMM, p-value = 0.77), suggesting the whales were feeding at the same rate throughout the study area. There was suggestive but inconclusive evidence that bouts were more dispersed for fin whales compared to blue whales after accounting for the distance to shore (GLMMzi p-value = 0.093), but there were no inter-species differences in the effect of distance to shore on both the distance to the closest bout (GLMMzi interaction p-value = 0.80) or the number of lunges per dive (GLMM interaction p-value = 0.96).

10.5 1,188 228 106 47 1,037

In two cases, a tagged whale passed through an area without stopping, where another tagged whale of the same species was feeding. Blue Whale # 2015\_4177 (a male) passed through an area 1 day after Whale # 2015\_840 (of unknown sex) had fed nearly continuously during daylight hours in the same location and more broadly for 12 d (**Supplementary Figure 6**). Similarly, fin

Whale # 2014\_5685 (a male) passed through a different area <1 day before Whale # 2014\_5838 (a female) spent 4 d feeding there (not shown). Both whales passing through were male and their movements were part of a larger circuit of southern California waters (**Figures 2**, **3**). The males' dive records were visually reviewed to ensure no feeding lunges had been missed by the lunge detection algorithm during this period of time.

One pair of blue whales each from 2014 to 2015 were recorded feeding in the same geographic space, in one case showing strong differences in dive behavior and in the other showing strong similarities. Whale # 2014\_5650 (a male) and Whale # 2014\_5803 (a female), tagged off Point Mugu, were in close proximity 13 times across an 8-d period, including times with both feeding and non-feeding behavior (**Supplementary Table 2**). While the characteristics of close-proximity periods were highly variable, Whale # 2014\_5803 generally dove more deeply and recorded more feeding lunges than Whale # 2014\_5650. On two occasions it dove to over twice the depth as Whale # 2014\_5650, and in one instance, was foraging when Whale # 2014\_5650 was not (**Supplementary Table 2**). During one of the longer periods of close proximity, both whales appeared to have been behaving similarly after passing closer than 0.5 km from each other (**Supplementary Figures 7, 8**), with both whales making foraging dives to a similar depth. An hour later, Whale # 2014\_5803 was feeding at almost twice the depth as Whale # 2014\_5650 (**Supplementary Figures 7, 8**) while in the same geographic space. In contrast, in 2015, Whale # 2015\_4177 (a male) and Whale # 2015\_5650 (also a male), tagged off San Miguel Island, fed in close proximity for 4.5 h on 8 July (the day they were tagged). These whales fed at approximately

the same depth (mean maximum dive depth = 273 and 212 m, respectively; **Supplementary Figure 9**) and shifted their feeding depth shallower during and after sunset with remarkable synchrony. Three other close-proximity events were recorded in 2015 but all were very limited in duration and number of dives recorded, so will not be presented. Other tagged whales from 2014 may have been in close proximity but the smaller number of Fastloc GPS locations collected by those tags did not

Whale # 2014\_5790 was not recovered and therefore no feeding data are available, but locations received through Argos are shown. The tracks for two additional fin whales of unknown sex are presented in Supplementary Figure 4.

allow for fine-scale resolution of their movements. In contrast to blue whales, there were no instances of close proximity between tagged fin whales, or between blue and fin whales, despite both being tagged in the same areas; in one case on the same day (e.g., blue Whale # 2014\_5655 and fin Whale # 2014\_5685), and in two cases of fin whales (Whale #s 2014\_5790 and 2014\_5838) being TABLE 2 | Summary of dives occurring during feeding bouts made by eight blue whales instrumented with ADB tags off southern California in August 2014 and July 2015.

*Total bouts* = *242.*

*Feeding bouts are sequences of dives with no more than 60 min between dives with recorded feeding lunges. Unknown sex whales are cases where no biopsy sample was collected. Note that scientific notation is used for the area of bout column, as values spanned several orders of magnitude.*

tagged on a day when a tagged blue whale (Whale # 2014\_5655) was re-sighted nearby.

#### DISCUSSION

The data presented here constitute the longest continuous dive behavior records collected to date for blue and fin whales, which allowed us to examine fine-scale behavior over timescales that previously have not been possible. Both species were observed and tagged concurrently during this study and subsequently occupied generally similar areas during their tracking periods, allowing for a comparative analysis between these two species, which occur in sympatry off California. Our focus was on characterizing feeding effort and its scales of variability over periods of multiple weeks. Short- to intermediate-duration (<6 h) feeding bouts were most numerous for both species and there was a positive correlation between bout duration and the number of feeding lunges made per dive within a bout. Blue whales have been shown to adjust their behavior and number of lunges made per dive based on the density of prey in the area (Goldbogen et al., 2015; Hazen et al., 2015), so the correlation between bout duration and number of lunges per dive suggests that the whales tracked in this study left lowerdensity prey patches quickly, while staying longer, and foraging more intensely, in higher-density patches. Some of the observed short-duration bouts may also represent the whales exhausting a highly localized abundance of prey. These results suggest the whales were following the marginal value theorem (Charnov, 1976), a foundational model of ecological theory that postulates that animals feeding in a patchy environment make decisions about when to depart a patch based on their assessment of its value, with the main prediction being that animals should spend more time in patches of higher quality (Mcnair, 1982).

Feeding bouts were more dispersed further offshore, although the whales were able to feed at the same rate throughout southern California waters, as the number of lunges per dive within a feeding bout did not change as a function of distance to shore. While there was evidence of inter-species differences in relation to distance from shore, they may have been the result of the more extensive use of offshore waters or more limited sample size of fin whales. The greater distance between bouts offshore suggests there may be two feeding strategies, with whales concentrating on highly localized, physically forced upwelling centers nearshore, such as off San Miguel Island (Fiedler et al., 1998), and more dispersed areas of elevated productivity offshore, TABLE 3 | Summary of dives occurring during feeding bouts made by three fin whales instrumented with ADB tags off southern California in August 2014 and July 2015.

*Total bouts* = *59.*

*Feeding bouts are sequences of dives with no more than 60 min between dives with feeding lunges. Unknown sex whales are cases where no biopsy sample was collected. Note that scientific notation is used for the area of bout column, as values spanned several orders of magnitude.*

likely driven by nearshore productivity that has been advected offshore, or by open-ocean Ekman pumping (Rykaczewski and Checkley, 2008). This has a range of implications for both habitat modeling and abundance estimation, as well as for managers trying to mitigate anthropogenic interactions, as the multiple scales of behavior and occurrence should be accounted for.

The spatial scale of feeding bouts was highly variable within and between individuals. The size of the feeding bout areas was likely an overestimate, as convex hulls are sensitive to irregular, concave shapes of the underlying points (Burgman and Fox, 2003; **Supplementary Figure 1**) and GPS locations were somewhat sparse in some cases, such that the number of locations may have been insufficient to define the true extent of the area being used for feeding. Despite this, median feeding bout size for fin and blue whales (5.6 and 7.6 km2 , respectively) appeared to correspond with the spatial scale of krill patches described off central California (1.8–7.4 km; Santora et al., 2011a). While some of the feeding bouts were significantly larger than the median values (>200 km2 ), krill patches up to 18 km in extent have been recorded in some years (Santora et al., 2011a), suggesting the larger feeding bouts may have been indicative of broad-scale areas of elevated krill abundance.

Logistics and expense limit direct study of key prey resources such as krill to fine-scales, from which broader-scale predictions are made and sometimes compared to the distribution of predatory species for validation (Santora et al., 2011a,b, 2014, 2017). Modeled krill distributions have also been used to better understand spatial and temporal aspects of their patchiness (Dorman et al., 2015; Messie and Chavez, 2017). The feeding data for tagged blue and fin whales presented here can offer further insight into the distribution and scale of prey patches across southern and central California waters, similar to how seabird foraging tracks have been used to infer prey availability and patch quality in other areas (Chimienti et al., 2017). Since direct observations of prey abundance were not available for this study, linking blue and fin whale feeding bouts to modeled krill distribution could be an additional step to both validate krill models and tie the overserved whale behavior more directly to prey abundance.

Tagged whale feeding behavior characteristics (like maximum dive depth or lunges per dive) fit into broadly similar ranges across individuals, although there was also evidence of fine-scale variability among individuals within those ranges, exemplified by the instances where two tagged blue whales (Whale # 2014\_5803 and Whale # 2014\_5650) were feeding in close proximity to one another but at different depths. Without knowing the structure of the prey field, it is difficult to be sure if these differences were related to the individual or the composition of prey being exploited. The observed differences may be a reflection of different individuals having different energetic requirements, allowing some whales to forage less intensively on lower prey concentrations (e.g., less dense prey at shallower depths), different age classes, or different prey species. However, across their entire tracks, Whale # 2014\_5803 fed at deeper depths than Whale # 2014\_5650, and this trend continued when the two occupied the same geographic space, suggesting the observed variability in dive behavior between individuals was likely due to different foraging strategies. It is not unusual for individuals to use a subset of their population's ecological niche due to variations of intrinsic traits or trade-offs that restrict an individual's ability to generalize (Bolnick et al., 2003). Further work is needed to better understand individual differences in rorqual behavior.

Previous work has found evidence of broad-scale spatial segregation between blue and fin whales off southern California (Fiedler et al., 1998; Irvine et al., 2014; Sirovic et al., 2015; Scales et al., 2017). However, our results showed that feeding behavior of tagged blue and fin whales was similar over a broad scale, with feeding occurring in relatively localized areas (median bout duration = 3.3 and 2.7 h, respectively; median bout area = 5.6 and 7.6 km2 , respectively) and with similar non-feeding periods in between (median time between bouts = 5.7 and 3.1 h, respectively). At finer scales, the characteristics of individual feeding bouts were also similar between the two species, with the

exception that fin whales recorded a higher maximum number of lunges per dive than blue whales (5.2 vs. 3.6, on average), consistent with previously observed species-specific differences in feeding rates (Friedlaender et al., 2015).

The observed similarities in feeding behaviors between the two species are likely a result of the highly specialized nature of lunge-feeding behavior (Goldbogen et al., 2006, 2011, 2017; Cade et al., 2016) and suggests the whales were feeding on krill given their highly stereotyped behavior (Cade et al., 2016). This further suggests the two species may potentially compete for a similar prey resource. Variations in the target prey species (Fossette et al., 2017), and even life stage of the same krill species (Santora et al., 2010), have also been shown to affect the spatial distribution of sympatric whale species. There was suggestive evidence that fin whale feeding bouts were more dispersed compared to blue whales, so it is possible that fine scale differences in the timing or location of feeding between the two species may be the underlying driver of the previously observed broad-scale differences in spatial distribution (Fiedler et al., 1998; Irvine et al., 2014; Sirovic et al., 2015; Scales et al., 2017). However, without an understanding of the underlying distribution and demographics of the local prey resources, it is difficult to attribute a reason to any inter-species differences in the spatial distribution of fine-scale feeding we observed.

Three tagged whales (two blue and one fin; all males) made clockwise circuits of southern California waters with only limited Irvine et al. Scales of Rorqual Feeding Behavior

feeding effort, despite passing through areas where other tagged whales were feeding. Blue whales have been shown to adjust their dive behavior based on the density of prey in the area (Goldbogen et al., 2015; Hazen et al., 2015), so it is possible they encountered prey in insufficient concentrations for them to feed upon. More broadly, males and females of central-place foragers are known to have different sex-based foraging strategies arising from different energetic requirements during reproduction and resulting in differential exploitation of the offshore and coastal environment (Breed et al., 2006; Austin et al., 2019). A similar process could be driving the observed differences between male blue and fin whales. Alternatively, blue whales are known to engage in social behavior while in southern California (Lomac-MacNair and Smultea, 2016). Visual inspection of the tag records revealed that many dives made by individuals of both species during the circuits of offshore waters had the characteristics of dives made while whales are vocalizing (shallow, flat-bottomed, minimal acceleration; see Calambokidis et al., 2007; Oleson et al., 2007; Stimpert et al., 2015). Additionally, it has been suggested that the better acoustic propagation properties of offshore waters (Sirovic et al., 2015) make them advantageous for reproductive calls (Lewis and Sirovic, 2018). Thus, we speculate these circuits in southern California waters by males of both species were related to reproduction, although we caution that this is based on a rather small number of individuals. Little is known about blue and fin whale breeding behavior, but our interpretation implies that courtship, or at least advertising and searching for a potential mate, may begin while on the feeding grounds.

The potential difference in movement and areas occupied by males and females has the added implication of possible sex-based differences in exposure to regional stressors like anthropogenic sound (Goldbogen et al., 2013; DeRuiter et al., 2017) or ship strikes (Berman-Kowalewski et al., 2010; Redfern et al., 2013), which are of particular concern in the Southern California Bight.  If further studies confirm our observations about sex–specific differences in habitat use (at least within southern California), managers may need to consider mitigation strategies that explicitly address differential impacts based on sex.

The movements of whales on the feeding grounds are generally assumed to be driven by the search for prey, such that fine-scale behavior can be directly linked to broad-scale ecological patterns like distribution. The results presented here support this idea, showing whales of both species quickly left presumably poor-quality prey patches, with the implication that variable prey density drives local whale movements, while the existence of multiple, persistent hotspots of krill aggregation along the west coast of the USA drive the broad-scale seasonal movement pattern (Abrahms et al., 2019; Palacios et al., 2019). However, sex also appeared to influence movement patterns and habitat use, at least within southern California waters. While this result is based on a relatively limited number of individuals, it constitutes a substantial refinement to our understanding of blue and fin whale behavior on the feeding grounds. Continuing to link fine-scale behavior to broader-scale movement is an important area of research, which was aided by the comparatively long data records generated by the tags used in this study. Developing a better understanding of the underlying mechanisms driving broad-scale behavior is critical to better assist the recovery of these endangered whales.

### DATA AVAILABILITY

The data used in this study are published as a Movebank Data Repository under a Creative Commons Zero license (Irvine et al., 2019).

### ETHICS STATEMENT

Tagging was conducted under the authorization of National Marine Fisheries Service Marine Mammal Protection Act/Endangered Species Act Research/Enhancement Permit No. 14856 and Oregon State University Institutional Animal Care and Use Committee Permit No. 4495.