Tagging and Tracking of Seabirds on the West Coast of Scotland

3. Methods

3.1 Fieldwork timing and objectives

The timing of fieldwork on Hirta, St Kilda and on Lunga, Treshnish Isles varied slightly across the three study years, mainly in response to improved understanding of the timing of breeding of the two species at the study sites, and the most appropriate time to carry out tagging operations. The timing of fieldwork and objectives of each campaign are shown in Table 1. All fieldwork was undertaken under all appropriate licences and permissions from NatureScot, BTO, the landowners and the RSPB HPAI exceptions panel.

Table 1: Timing and objectives of fieldwork conducted on St Kilda and Treshnish Isles 2021-2024.

Leach’s Storm-petrel

Hirta, St Kilda 57.81⁰ N -8.60⁰ W

European Storm-petrel

Lunga, Treshnish Isles 56.49⁰ N, -6.42⁰ W

3.2 GPS tracking

The small size of Leach’s Storm-petrels and European Storm-petrels currently precludes the use of tags which allow remote data download. Some manufacturers are in the process of developing miniaturised GPS tags capable of remote download and a small number of prototype tags supplied by TechnoSmart were trialled during the 2023 fieldwork for this project. However, they performed very poorly and were not used further. Only one of four supplied tags collected data during pre-deployment trials and the fourth tag failed to download data to the base station during nine days of deployment. The bird was recaptured and the tag was removed. The tag antenna was found to have broken and the solar panel had suffered corrosion. Therefore, throughout this project we used archival GPS tags supplied by PathTrack (nanoFIX® GEO-Mini) which are the lightest (c. 0.93g) and smallest (18 x 11 x 4mm) commercially available avian GPS tags, and which had been used successfully for tracking European Storm-petrels previously (Bolton 2021). To maximise the chances of recapturing birds for tag retrieval and data download, individuals were tagged at sites where an active breeding attempt was confirmed or considered very likely (see below). Nests were found through playback of recordings of their calls from a handheld loudspeaker (Gilbert et al. 1998). This method was used both during the day and at night. Responses from birds in burrows indicated possible breeding sites, which were then monitored for potential use as study nests. Where possible, the status of the nest (e.g., pre-breeding, incubation) was confirmed through direct visual inspection of the nest chamber, or by use of an endoscope where the nest chamber was not directly visible. If the nest chamber could not be inspected, consistent occupation of a nest site indicated through repeated response to playback over multiple days was taken as an indication of breeding.

Where nest chambers were sufficiently accessible, birds were caught for tag deployment by direct removal from the nest by hand. However, the nest entrances of Leach’s Storm-petrels were too narrow and long for removal by hand, so birds were caught by placing a trap at the nest entrance to catch birds as they departed to begin a foraging trip. GPS tags were attached to the central four tail feathers using thin strips of adhesive Tesa® tape (c. 0.1g). The tag attachment method was temporary and if it was not possible to recapture the bird to remove the tag, it would fall off after a period of a few weeks. To retrieve tags deployed on Leach’s Storm-petrels, we placed a trap at the nest entrance to catch the tagged bird on its return. Traps were constantly observed using thermal imaging equipment (Helion Pulsar binoculars or monocular, various models) for capture of returning birds. For European Storm-petrels, where nest entrances were generally unsuitable for use of such traps, but nest chambers were accessible by hand, we checked nests daily following tagging and if the tagged bird had returned from a foraging trip and was present in the nest during the day, the tag was retrieved the following evening. During the post-brooding stage, when the tagged birds did not remain in the nest during the day, we watched nests using thermal imaging equipment at night until the tagged bird was observed to return, and the bird was caught by hand in the nest chamber to retrieve the tag.

Assessment of potential impacts of GPS tagging and tracking on (i) foraging behaviour and (ii) nest survival rates was complicated by the lack of accessibility of many nests chambers for marking birds (for individual recognition to assess nest attendance patterns) and assessment of nest contents (to assess egg and chick survival rates). Most Leach’s Storm-petrel nests could not be accessed directly to mark birds or view the nest occupants. However, the majority of nesting Leach’s Storm-petrels responded consistently to late-evening playback at the burrow entrance of recordings of calls. Since calls differ between the sexes, playback therefore permits identification of the bird present at each nest each day and allows collection of data on nest attendance patterns of untagged birds. In 2021 we compared foraging trip duration of tagged birds during the period of tag deployment with the foraging trips of: (i) the same individuals prior to tagging and (ii) all untagged birds. The first analysis was conducted using Wilcoxon matched-pairs analysis, and the second analysis used a linear mixed effects model. For the mixed model, we included all first trips of tagged birds (N = 14) as well as all trip duration data recorded at a further 37 nests (N = 94) throughout the study period (18^th June – 10^th July). As we did not have information on breeding stage for all nests, we could not include this as a random effect in the model, and instead included the date of the start of the trip as a fixed effect (Julian date, centred and rescaled). Nest was included as a random effect to account for multiple measurements from some nests. Trip duration data was square-root transformed to meet the assumption of normality of residuals. For both analyses we inferred foraging trip duration for untagged birds from the duration of periods between successive incubation bouts. Tracking indicated that birds spent the entire period between successive incubation bouts at sea.

For European Storm-petrels, the inaccessibility of many nest sites similarly meant that it was not always possible to individually mark birds or identify which of the breeding pair was present. Nest attendance data were therefore not sufficiently complete to identify differences in nest-attendance patterns between tagged and non-tagged birds before and after tagging. In 2021 daily nest survival was compared for the 23 nests where tagging of birds occurred and 52 control nests to check for an effect of tagging on nest survival. Whilst periods of temporary egg neglect are common in this species (Davis 1957), a period of six consecutive days of an egg being unattended was considered to indicate nest failure. The loss or damage of an egg or chick was also considered to be a nest failure. Nests were monitored for up to 27 days and we compared daily survival rates of nest where tagging occurred and control nests, using a binomial model.

3.3 GLS tracking

In order to track birds throughout the course of the non-breeding season, geolocator (GLS) tags (0.43g W30A9-SEA Migrate Technology) were deployed in July and August 2021 using long-term attachment methods, based on protocols developed by researchers in Spain (European Storm-petrel) and Canada (Leach’s Storm-petrel). Each tag was attached using monofilament cord and marine epoxy to an overlapped darvic ring, and fitted to the tibia. To minimise handling of birds and associated disturbance in the colony, most GLS tags were fitted during the same capture event as GPS tag retrieval. Additional birds were tagged with GLS devices after removal from nests by hand, either during the chick-rearing stage or late in the incubation stage, when birds are more robust to disturbance. Where possible, the mates of all tagged birds were ringed and processed as controls.

To retrieve GLS tags deployed in 2021 (and to assess any impacts of GLS tagging) all nests that had been occupied by birds fitted with GLS tags, and their controls, were monitored regularly when fieldwork resumed in 2023 using playback and/or visual inspection of the nest contents. Traps and/or mist nets were deployed at night at the entrances of nests that were found to be occupied, and also (on a less frequent basis) at those nests were no response to playback was ever obtained (in case occupants did not respond to playback). We noted the extent of recapture effort directed at both GPS-tagged birds and their controls, in order to avoid biasing effort towards either group. Any recaptured birds that have been fitted with GLS tags in 2021 (and controls) were inspected for sign of injury and any GLS tags removed. Potential impacts of GLS tagging were assessed by comparing the recapture rates of GLS tagged birds and controls.

3.4 Assessment of HPAI impacts

In 2023 blood samples were taken from the brachial (wing) vein 12 Leach’s Storm-petrels on St Kilda (Figure 2) and four European Storm-petrels on Treshnish, under a Home Office Licence issued to Mark Bolton, working in collaboration with Dr Emma Cunningham (Edinburgh University). Additional work was undertaken by RSPB (under separate funding) to resurvey the main breeding colonies of Leach’s Storm-petrels on St Kilda (Dun, and Carn Mor, which jointly held nearly 80% of the entire population of St Kilda in the most recent survey in 2018). This work was carried out as part of the RSPB HPAI seabird survey project. Data analysis is scheduled for 2025 and will be reported separately.

3.5 GPS Data Analysis

I. Data processing and trip metrics

Data were processed in RStudio (R Core Team 2023 & Posit team 2023). To provide information about trip metrics, GPS tracks were split into trips based on departure from and return to the colony using the ‘tripSplit’ function of the track2KBA R package (Beal et al. 2021). Departure and return buffer thresholds were set on the basis of visual inspection of tracks. For Leach’s Storm-petrels, a departure buffer of 5 km and a return buffer of 30 km was used. For European Storm-petrels a departure buffer of 5 km and a return buffer of 10 km was used because the tracks included movements of birds foraging relatively close to the colony. To validate the trips identified by the ‘tripSplit’ function, trips were also manually labelled in ArcPro. Care was needed to ensure ‘tripSplit’ correctly identified partial trips. An additional label of ‘Almost Complete’ was added for trips where, due to tag failure or low battery, the track terminated just short of the colony or there was a longer interval between fixes than the tag was programmed to record. These “Almost Complete” trips were included in summaries of maximum distance from the colony (where maximum distance was known) but not for total distance or trip duration (which were not accurately recorded).

We used the ‘tripSummary’ function of track2KBA to estimate the total distance travelled, maximum distance from the colony and duration for each trip. The distance between the first and last point of the trip and the colony location was added to distance calculations in instances where trips started or ended at sea. However, when estimating trip duration, the ‘tripSummary’ function assumes the start and end times of a trip are those recorded at the first and last GPS location for that trip. Where these locations are at sea (i.e. because the bird’s presence at the colony before and after the trip was not recorded, due to difficulty for the tag in viewing satellites when the bird is below ground), the trip duration will be underestimated. To account for this, tracks were rooted to the colony by appending the colony location at the start of the trip, one timestep prior to the first recorded location at sea. This ensured that trip duration and total distance travelled were comparable between trips. No additional timestamps were added to the end of trips to avoid overestimating trip durations. For Leach’s Storm-petrels, trips were labelled as being from one of two breeding stages: incubation and chick-rearing. Where there was evidence that disturbance during trapping attempts or a failed breeding attempt had altered the birds’ behaviour, trips were excluded from the overall summaries and further analyses, as follows. One bird lost its egg whilst at sea on a foraging trip, and thereafter did not remain at the nest by day. Instead it made a series of trips that included several short-range, short-duration trips, unlike other birds with active nests (Table 4). In addition, there were a few occasions when birds returned to the colony during the egg stage, but did not enter the nest to resume incubation – instead they returned to the sea for a further day, before returning to the nest the following night. This resulted in a few short duration (one day) and short range trips being recorded during incubation, (when birds are typically at sea for 2-3 days). These “aborted” nest visits occurred on nights when captures of these birds were attempted, Whilst it is not clear to what extent the observed behaviour was normal, or whether birds were hindered from entering their nests by nocturnal fieldwork activities, as a precaution these short trips following are excluded from analysis. These cases applied to the following trips for Leach’s Storm-petrels during the incubation phase: 2021_BT27403_44_02 (Nest Failure), 2021_BT27403_44_03 (Nest Failure), 2021_BT27403_44_04 (Nest Failure), 2021_BT27403_44_05 (Nest Failure), 2021_BT27403_44_06 (Nest Failure), 2021_BT27403_44_07 (Nest Failure), 2021_BT27406_38_02, 2021_BT27420_46_02, 2023_BT27462_45_02, 2024_NT63501_122_02, 2024_NT63501_122_03.

To assess the degree of individual specialisation in space use we used the function ‘indEffectTest’ from track2KBA (Beal et al. 2021). The function uses a bootstrapping procedure to compare random pairs of trips which comprise either (i) two trips from the same individual, or (ii) two trips from different individuals. An asymptotic two-sample Kolmogorov-Smirnov test compares within-individual and among-individual variance in 50% UDs. To allow for the possibility that foraging areas might differ between breeding stages, the test was applied separately to foraging locations during incubation and chick-rearing stages (i.e. pairs of trips from different individuals were drawn from the same breeding stage).

For European Storm-petrels trips from 2021 only were summarised by breeding stage as all birds tracked in 2023 and 2024 was conducted during the chick-rearing period. Where we found evidence of individual specialisation in space use, such that multiple trips performed by the same individual could not be regarded as independent, or broadly representative of a large sample of tracked individuals, summary statistics were derived from the output of linear mixed-effects models with individual ID used as a random term to derive individual-level (not trip-level) estimates. Since model outputs relate to standard error (SE) of the mean rather than sample standard deviation (SD) we estimated sample SD from equation 1 below (DF = Degrees of Freedom).

SE x √(DF + 1) - eqn 1.

Where there was no evidence of individual specialisation, we considered that multiple trips by the same individual were representative of a sample of single trips from multiple individuals, and calculated summary statistics (and conducted statistical hypothesis testing) directly on the raw trip-level data. Where there was evidence of individual specialisation we calculated summary statistics and conducted hypothesis testing using a random effects model to account for non-independence of trips undertaken by the same individual.

The ‘lme4’ package (Bates et al., 2015) was used to build linear mixed-effects models with individual as a random effect. For each variable (trip distance, maximum distance and duration) four models were produced. The null model contained only individual ID as a random effect, this was compared with a model including both year and breeding stage, just year or just breeding stage. Where multiple models revealed a significant effect, AIC was used to select the best performing model.

Representativeness was assessed using the ‘repAssess’ function in the track2KBA R package (Beal et al. 2021). This function assesses how representative is the assessment of space use derived from the tracked sample, compared with the estimated extent of colony-level space use. The function implements a bootstrapping procedure to iteratively resample the utilisation distributions of individuals within the sample and then models the relationship between space use of all individuals in the sample and sample size. The output is a score between 0 and 100 reflecting the extent of colony-level space use identified from the tracked sample. Scores above 70% indicate that the dataset is highly representative of the space use at the colony level (Beal et al. 2021).

Maps of bird movements were produced in R using the ggOceanMaps package (Vihtakari 2024). Manual labelling of trips to validate the automated labelling of trips by track2KBA and production of the maps showing overlap between the European Storm-petrels and OWF lease areas was performed in Arc GIS pro 2.8.8 (ESRI 2021). Data on the location of Scottish Offshore Windfarm Areas (OWFs) and those in English, Welsh and Northern Irish waters are freely available via the Crown Estate spatial hub hosted on Arc GIS online (Crown Estate 2021). Data on the location of relevant OWFs in Republic of Ireland waters were sourced from Marine Renewable Energy Ireland (MaREI 2021) via the Irish government’s open data portal (Open Data Portal 2024).

II. Home Range Estimation

The observed utilisation distributions (UDs) of both species were estimated using the ‘ks’ package in R to produce 25%, 50% and 95% kernels (Chacon & Duong (2018). UDs were produced across all birds, by year, by time of day and by breeding stage. The ‘ks’ package has the benefit of including a number of ‘bandwidth selector’ functions for automating the selection of the optimal smoothing factor for building the kernel density surface. For our data, the ‘Hpi’ function was selected as the bandwidth selector. This was chosen because it can handle the inclusion of negative values associated with co-ordinates in decimal degrees.

Marine Habitat Association

To help inform future spatial planning of marine development such as offshore wind, a key objective of this project was to improve our understanding of marine habitat use by foraging Leach’s Storm-petrels and European Storm-petrels. A common approach to understanding habitat use is resource selection modelling whereby, for each observed location, a number of pseudo-absences are generated within a given domain (Johnson 1980, Boyce et al. 1999, Aarts et al. 2012). The domain refers to the area available to the animal, i.e. where the animal could feasibly have been present, and is set by the researcher. Resource selection functions (RSF) then compare the environmental conditions in these available locations with the conditions at the locations where the animal was observed to be present (Paton and Matthiopoulos 2016). In the context of the standard resource selection function there is no single optimal method for determining the area available to an animal (Paton and Matthiopoulos 2016). Moreover, the exact choice of the available domain dictates the scale at which habitat selection is being considered and is therefore key to model interpretation (Johnson 1980, Northrup et al. 2022). For seabirds, common approaches include using the maximum observed distance from the colony (e.g. Wakefield et. al. 2017) or the observed home range (Cleasby et al. 2024). These methods have been shown to effectively characterize broad-scale seabird habitat use, providing reasonable predictive performance (Wakefield et al. 2017). However, RSF approaches can introduce issues arising from the non-independence and spatio-temporal autocorrelation of location data obtained by telemetry because they are derived from methods developed for point count data (Alston et al. 2023). To address some of these issues and make use of the fine-scale movement data from telemetry studies, Fortin et al. (2005) suggested resampling step lengths (distances between successive observed locations) and turn angles (deviations from previous bearings) to generate random movements and thus available points, conditional on the previously observed locations. This process creates stratified datasets, with a distinct set of available points associated with each observed location. Such models are commonly termed step selection functions (SSF) and are focussed on fine-scale habitat selection at the scale of GPS sampling intervals (Klappstein et al. 2021, Alston et al. 2023).

Here, we implemented a SSF procedure to investigate fine-scale habitat selection in European and Leach’s Storm-petrel (Mercker et al. 2021). Step selection functions are well suited to the analysis of habitat selection in wide ranging species such as petrels. This is because SSFs consider an animal’s track as a series of discrete sequential steps separated by constant time intervals, where each step is defined by the distance moved and turning angle (Thurfjell et al. 2014). The habitat characteristics associated with landing point for the next step are compared with those from a set of pseudo-absences generated from the observed range of turning angles and step lengths as illustrated in Figure 3. A key advantage of the SSF approach is that it resolves the challenge associated with defining what is ‘available’ to an animal by using parameters of animal movement to define availability (Forester et al. 2010). An SSF approach enables effective characterisation both long-range attraction and shorter-range local resource selection.

Collating Suitable Environmental Covariates for SSF

To examine habitat selection in both storm-petrel species, we chose to incorporate a range of potentially important environmental covariates to which storm-petrels might respond. Given the wide array of covariates that could be considered, we focused on those that have been shown to be important for other UK seabirds (e.g. Wakefield et al. 2017) to generate a suitable yet manageable set of environmental covariates for modelling purposes. Additionally, covariates had to be freely available over the same spatio-temporal extent as the storm-petrel tracking data, which further reduced the number of available covariates.

In total, the static environmental covariates we selected were: (1) Distance by sea from the colony (km), where distance by sea was calculated using the centre point of each cell contained within a regular grid (1 km² resolution) encompassing the area defined as available to each colony. This area was defined as a radius around the colony that extended to 1.1 x the maximum foraging range observed at each colony (Wakefield et al. 2017); (2) water depth (metres, m) which was derived from Global Bathymetry and Topography at 15 Arc Sec: SRTM15+ V2.5.5 bathymetry data (Tozer et al. 2019) with a nominal resolution of 15 arc seconds (approximately 500 x 500 m pixel size at the equator); (3) seabed slope (°) which was generated from the same bathymetry data used to estimate water depth using the terra package; (4) Distance from the coast (km) at 1 km x 1 km resolution; (5) benthic habitat derived from the EUSeaMap 2021 (Vasquez et al. 2021). Data from EUSeaMap was processed in R to simplify the 29 habitat categories into six broad categories: ‘1:Rock’, ‘2:Coarse Substrate’, ‘3:Sediment’, ‘4:Sand’, ‘5:Mud’, ‘6:Reef or other biogenic habitat’ and ‘Unlabelled’. The numeric code (1-6) associated with each habitat category was then used to generate a raster surface at a resolution of 0.01° x 0.01° (Approximately 1.11 km x 1.11 km). In addition, the following dynamic variables included were: (6) Sea Surface Temperature (SST, °C) at 1 km² resolution using Multi-Scale Ultra High Resolution (MUR) SST data (MUR SST was accessed on 15/02/2025). Daily sea surface temperature (SST) data were downloaded for each date on which tracking data were available across the years. For each spatial grid cell, the average SST for each year of the study was calculated using the daily SST data corresponding to the specific tracking period in that year (from the date of the first to the last GPS observation) ; (7) the magnitude of ocean SST front gradients (°C / 1 km) estimated via the grec package (Lau-Medrano 2020) using daily SST data and the front detection algorithms developed in Belkin and O’Reilly (2009). As with SST, yearly averages of thermal front gradient were calculated by averaging across daily values that covered the tracking period in each year; (8) mass concentration of chlorophyll-a in sea water (CHL-a, mg/m³) at ~7 km² resolution monthly from June to August, downloaded from the Copernicus Marine Service (CMS)(CMS, date accessed: 22/01/2025). Monthly maps of chlorophyll each year were then averaged to give an estimate of mean chlorophyll per spatial grid cell across the summer (June to August) in that year; (9) the magnitude of ocean chlorophyll gradients (mg/m³ / 7 km) estimated via the grec package (Lau-Medrano 2020) using monthly chlorophyll data and the front detection algorithms developed in Belkin and O’Reilly (2009); (10) Simpson-Hunter stratification parameter (S, Simpson et al. 1974) this index (S = log₁₀(hu⁻³)) combines measurements of u = mean current speeds (ms^-1) and h = depth (m) into a single measurement which helps to identify mixed, frontal, and stratified water columns (Bowers & Simpson 1987, Holt & Umlauf 2008). Water depth data was obtained as described above and current speed data was obtained from the Copernicus Marine Service (CMS, date accessed: 07/02/2025). Current speed was downloaded for each date for which tracking data was available at 12:00:00 (midday) at ~ 8 km² resolution and averaged across the observed tracking period each year. To calculate the Simpson-Hunter parameter, the current data were smoothed using bilinear interpolation to ensure its resolution matched that of the depth data; (11) Net primary production of biomass expressed as carbon per unit volume in sea water (NPP, mg/m³/day) at a monthly ~ 7 km² resolution (CMS). For each year, average NPP per grid cell over the tracking period was then calculated.

Fitting SSFs and Model Selection

For both Leach’s and European Storm-petrel, step-selection functions were fitted using the hmmSSF R package that allows for the inclusion of non-linear habitat associations (Klappstein et al. 2021). One requirement of SSFs is that movement data are sampled at regular time intervals to ensure step lengths and movement remain consistent throughout a tracking period and across different individuals. Movement trajectories were regularised for each identified foraging trip in turn to ensure periods at the colony in between trips. To regularize movement trajectories, we used the track_resample function from the amt package to resample tracks at a regular 1 hour interval with a tolerance of 1 minute (min) for both petrel species. We chose a 1-hour interval as this reflected the lowest sampling resolution observed for some individuals in both datasets. Some individuals were tracked at higher resolutions than 1-hour intervals (15 min or 30 min intervals), however resampling at such shorter sampling intervals would lead to loss of information from individuals with coarser sampling regimes or an increasing reliance on interpolated points for such individuals.

For each observed point, we generated a paired sample of 50 potentially available points. These available points were sampled using a non-parametric movement kernel within the hmmSSF package. Specifically, available points were selected uniformly from a disc centred on the previous observed location, with the disc radius (R) set close to the maximum observed step length to ensure adequate spatial coverage (Klappstein et al. 2022; Michelot et al. 2024). The advantage of using a non-parametric kernel to generate available points is that we do not need to provide a specific distribution for either step lengths or turning angles. However, the non-parametric approach requires more available points to be sampled for accurate parameter estimation hence our choice to sample n = 50 available points per observed point (Klappstein et al. 2024). Once a dataset of observed and available points was constructed, we then used the terra package to extract the relevant environmental covariates observed at each point.

SSFs were fitted as conditional logistic regression models implemented as generalised additive models (GAM), which allow covariates to be fitted as non-linear smoothers, using the mgcv R package. Following, Klappstein et al. (2024) these GAMs were constructed as Cox Proportional Hazard (CPH) models. Due to the relatively large number of environmental covariates examined we conducted model selection in a forward stepwise fashion. First, we ran models including each candidate environmental covariate in turn as either a standard linear effect or a non-linear smoother and judged how well model fit was improved compared to a null model that contained no environmental covariates. For some variables which took large values on their original scale (distance from colony, depth, NPP, distance from the coast) we also considered models in which model coefficients and non-linear smooths were estimated on the log-scale. Following Klappstein et al. (2024) our null model still included terms for step lengths and turning angles modelled using non-linear splines to describe the distribution of these movement variables. All smoothers were fitted as shrinkage splines (so that smooths can be penalised out of the model if not required), except for those related to turning angles, where cyclical splines were used to account for the circular nature of this parameter. Due to potential collinearity between covariates, we excluded models that included two covariates with a correlation greater than r = 0.4. In such cases, the covariates selected first in the stepwise modelling process were preferred. Additionally, we did not include both water depth and the Simpson-Hunter parameter in the same models, as water depth is used as the denominator in the equation to estimate the latter. In the results we report the finally selected model for each species.

Model performance was assessed using model Akaike Information Criterion (AIC) scores, by default the mgcv package reports a ‘corrected’ conditional AIC which accounts for uncertainty in model smoothing parameters which will tend to favour simpler models (Wood et al. 2016, Klappstein et al. 2024). The best performing model in this first step was then selected as the basis for the next step in the model selection process until reaching a model at which including further covariates no longer led to improvements in model AIC. Once the best performing model was selected, we considered the effect of including additional two-way interactions between the environmental covariates included in this model and the following group-level predictors: 1) whether an observation was recorded at day or night (calculated using the r suntools package (Bivand & Luque 2024) to determine the timing of sunrise and sunset at specific locations and timestamps); 2) whether a bird was classified as incubating or chick rearing; and 3) the year of the study. If the inclusion of a two-way interaction improved model fit, it was retained otherwise it was removed.

3.6 Geolocator data processing

Geolocation data processing was carried out using R 4.2.3 (R Core Team 2023). We used the preprocessLight function of the ‘TwGeos’ package (Lisovski et al. 2016) to process log-transformed light data using a light threshold of 1 log lux. This threshold allowed us to define sunrise and sunset (twilight events), which occur when the light values exceed or fall below this threshold, respectively. With this function, we estimated the hour of sunrise and sunset, inspected the integrity of the light curve of each day and manually adjusted the time of sunrise or sunset where the time appeared incorrect in comparison with the observed light profile of preceding and subsequent days. These incorrect times likely occurred due to interference of feathers covering the light sensor. We manually adjusted the incorrect time of twilight to the mean time between the previous and the following days, except in those cases when there was a small peak of light (below the light threshold we defined) before (or after, in the case of sunset) the twilight time assigned automatically; in those cases, we moved the transition to this small peak of light. Subsequently, we used the ‘Solar/Satellite Geolocation for Animal Tracking (SGAT)’ package (Sumner et al. 2009, Lisovski & Hahn 2012), which applies Markov chain Monte Carlo (MCMC) simulations to estimate and refine the locations of the tracked individual. For this analysis, we first calculated the zenith angle and the error distribution around the twilights based on the calibration period of the geolocator at a known location. We also generated a gamma distribution of flight speeds between 0 and 25 km/h based on the average speed values from the GPS tracking element of this project. We then used the thresholdPath function to obtain the initial path of the bird, which is needed to begin the MCMC simulations. Furthermore, we generated a spatial mask to avoid any locations over land. We used these parameters and started by drawing an initial 200 samples for burn-in and tuning of the proposal distribution. A further 300 samples were drawn to evaluate chain convergence before drawing another four MCMC chains of 3,000 samples each to describe the ultimate posterior distribution and the most likely migration path. Despite using a high tol (tolerance on the sine of the solar declination) value when running the thresholdPath function (which defines how many locations should be linearly interpolated around the equinox), we still had to remove unrealistic values near the equinox periods.

3.7 Overlap with Offshore Wind Areas

For European Storm-petrels, a time in area analysis was performed to estimate the proportion of time spent during each foraging trip in proposed OWF areas (no foraging trips of Leach’s Storm-petrels overlapped with OWF lease areas). To prevent underestimation of the number of individuals interacting with OWFs and control for different sampling frequencies, tracking data were linearly interpolated to a 1-minute sampling interval. This was achieved using the ‘mt_interpolate’ function from the move2 package in R (Kranstauber et al. 2024).

Spatial data on the locations of OWFs in British and Irish waters was used to generate a raster surface of unique IDs whereby each ID related to a specific wind farm in the data set. The ‘terra’ package (Hijmans 2025) was then used to overlay the tracking data (including the points generated from interpolation) onto the raster surface to extract the ID of any OWFs. Points which did not overlap with an OWF were assigned a value of zero. The numeric ID was then used to append the name of each relevant OWF to the data set. The total time spent in each OWF per trip was then used to calculate the percentage of time spent in OWFs during each foraging trip.