The Use of Animal-Borne Biologging and Telemetry Data to Quantify Spatial Overlap of Wildlife with Marine Renewables

Abstract

The growth of the marine renewable energy sector requires the potential effects on marine wildlife to be considered carefully. For this purpose, utilization distributions derived from animal-borne biologging and telemetry data provide accurate information on individual space use. The degree of spatial overlap between potentially vulnerable wildlife such as seabirds and development areas can subsequently be quantified and incorporated into impact assessments and siting decisions. While rich in information, processing and analyses of animal-borne tracking data are often not trivial. There is therefore a need for straightforward and reproducible workflows for this technique to be useful to marine renewables stakeholders. The aim of this study was to develop an analysis workflow to extract utilization distributions from animal-borne biologging and telemetry data explicitly for use in assessment of animal spatial overlap with marine renewable energy development areas. We applied the method to European shags (Phalacrocorax aristotelis) in relation to tidal stream turbines. While shag occurrence in the tidal development area was high (99.4%), there was no overlap (0.14%) with the smaller tidal lease sites within the development area. The method can be applied to any animal-borne bio-tracking datasets and is relevant to stakeholders aiming to quantify environmental effects of marine renewables.

1. Introduction

Marine renewables (fixed and floating offshore wind; wave; and tidal energy) have been identified as key to the transition to net-zero emissions [1,2]. In order for consenting decisions for marine renewable energy developments to be made responsibly, any potential negative impacts on key ecological receptors (birds, fish, marine mammals, turtles) need to be assessed [3,4,5]. A key constraint for marine renewables is understanding how animals use development sites, and the risk of collision that they may be exposed to [5,6,7,8]. Taking tidal energy as an example, fatal collision with dynamic parts of operating devices (e.g., rotating blades) has been identified as a major concern for marine vertebrates [5,7]. Due to the small number of operational tidal turbines in the water worldwide, detecting collision events, or the absence thereof, remains challenging [3]. Where such information is lacking, an alternative is to quantify space use of animals in tidal stream environments [9,10,11]. Assessing spatial overlap of these areas with areas in which turbines are likely to occur, that is tidal energy development sites, can then inform about the level of risk to wildlife [5,12,13]. Any connectivity with renewable energy developments will then not only identify populations at risk of exposure but also inform the need for and extent of any additional research. This is especially relevant as the marine renewable energy sector continues to expand, especially offshore wind [14,15].

Animal-borne biologging and telemetry devices are used extensively in marine vertebrate research to record movements above, at, and below the sea surface [16,17,18]. Telemetry devices such as GPS loggers can record the location (typically to within <10 m) of an individual animal down to every second, thereby providing accurate data at fine spatiotemporal scales [19,20]. Time-depth recorders (TDRs) are examples of biologging devices that continuously record underwater pressure from which important dive metrics (e.g., profile, depth, duration) can be extracted [21,22]. Such data have previously been used in marine conservation and management, for example to help designate and evaluate effectiveness of Marine Protected Areas or to identify interactions between fisheries and vessels [23,24,25,26,27,28]. At a finer scale, simultaneous deployments of biologging and telemetry devices on mobile receptors have the potential to identify areas of spatial overlap with marine renewable energy developments [12,29,30,31].

The extraction of utilization distributions from animal-borne tracking data to identify areas used by animals is a common method in ecology, e.g., [32,33,34,35,36]. The outputs are intensity of area use polygons that can be readily interpreted and are highly comparable between sites and species [37]. Advantageously, biologging devices supplement GPS positions with metrics relevant to marine renewables such as flight height or dive depth [38,39]. This allows for the stressor of interest (e.g., collision with rotating blades) to be explicitly considered in the generation of utilization distributions.

Despite the potential of animal-borne tracking data, a number of factors inhibit their use in the assessment of effects from marine renewables. These include ethical and logistical considerations inherent to deployment of tags on wild animals as well as complexity in the processing and analysis of these types of data [17,40,41]. In order for this high-resolution data to be useful to stakeholders within the marine renewable energy sector, there is therefore a need for straightforward and reproducible workflows, including the development and presentation of relevant code [42]. Consequently, the aims of this study were: (1) to develop an analysis workflow in order to rapidly process high-resolution biologging and telemetry tracking (point) data and produce areas of use (polygon) data, for incorporation in environmental impact assessments or for decision-making purposes regarding future siting of marine renewable energy developments; and (2) to demonstrate the methodology for European shags (Phalacrocorax aristotelis), one of the species assessed to be most vulnerable to potential adverse effects, collision in particular, from tidal turbines [8,43], including quantifying spatial overlap with a major tidal energy development area.

2. Materials and Methods

2.1. Study Site

The Pentland Firth is a channel between the north coast of mainland Scotland and the Orkney archipelago, connecting the North Atlantic Ocean and the North Sea (Figure 1). The extremely fast tidal currents in this channel (>5 ms⁻¹ at mean spring tide) have led to it being an important area for tidal development in the UK and to its designation as a tidal draft plan option [44,45,46] (https://www.renewableuk.com/page/WaveTidalEnergy, accessed on 1 July 2020). Within the Pentland Firth, there are two tidal lease sites, the Inner Sound and Ness of Duncansby. The Inner Sound is currently leased by MeyGen Ltd. and since 2016/2017 contains four grid-connected horizontal axis tidal turbines [47]. Spatial polygons of tidal draft plan options and tidal lease sites were obtained from Marine Scotland Science (http://marine.gov.scot/information/draft-sectoral-marine-plans-wind-wave-and-tidal-2013, accessed on 1 July 2020 [14]) and Crown Estate Scotland (https://www.crownestatescotland.com/, accessed on 1 February 2020).

2.2. Bird-Borne Biologging and Telemetry Dataset

Data from European shags simultaneously carrying GPS tags and TDRs were obtained through the Environmental Research Institute and RSPB FAME/STAR projects (www.fameproject.eu, accessed on 1 October 2018). The dataset consisted of six birds tagged within the Pentland Firth during the breeding season (May–June; Five on Muckle Skerry (2012 n = 2; 2013 n = 1; 2014 n = 2) and one on Stroma (2012 n = 1). (Figure 1,

Table 1. Summary statistics of the six shags included in the analysis.

ID	Start Date	End Date	# of Fixes	# of Trips	Duration (h)	Mean ± SD Distance to Colony (km)	Max Distance to Colony (km)
EUSH616	2012-05-26	2012-05-26	236	1	10.8	11.1 ± 4.2	16.1
ORK0445	2012-05-21	2012-05-24	328	5	59.0	4.5 ± 2.7	8.4
ORK0451	2012-05-20	2012-05-23	429	7	71.0	1.2 ± 0.5	3.1
ORK0694	2013-06-14	2013-06-17	184	6	65.2	0.6 ± 0.1	0.9
ORK0795	2014-06-09	2014-06-11	127	3	43.6	1.8 ± 0.5	3.2
ORK0797	2014-06-09	2014-06-10	80	2	24.3	2.3 ± 0.7	3.4
Total			1384	24	273.9	3.6 ± 1.45	16.1

, Appendix A). The combined weight of devices was always <3% of shag body weight, and further details on bird handling, attachment method and devices can be found in [48,49,50].

2.3. Workflow Description

2.3.1. Preparation and Processing

Data preparation, processing and analysis were carried out in R version 3.6.1 and R Studio version 1.2.1335 [51]. The steps for preparation and processing, adapted from [27], are outlined in Figure 2 and followed methods outlined in [48]. The simultaneous GPS and TDR datasets from each individual were merged to the nearest GPS point in time, resulting in a dataset with dive locations identified. The dataset was then rediscretized (linearly interpolated) to 120 s intervals between fixes using the adehabitatLT package [52].

During the breeding season, shags are a central place foragers and are therefore constrained to a maximum distance from the nest, making trips from the nest to and from foraging sites in order to defend, incubate and feed their offspring [53]. Foraging trip classification is therefore an inherent part of breeding seabird track analysis. Foraging trips were defined as a minimum 30 min round trip from and to a 300 m buffer around the nest. These thresholds were chosen to eliminate resting and washing behaviors as far as possible [54,55]. As the stressor of interest was collision with tidal energy devices which are deployed at, in the case of floating tidal, or below the sea surface, GPS fixes likely to be flight as well as those on land were removed, leaving only at-sea fixes likely to be part of foraging (dives or sitting on the sea surface). Flight fixes were identified based on a speed and tortuosity threshold of groundspeeds over 5 m s⁻¹ and turning angles over 45 degrees [55,56]. The speed threshold was chosen to account for the exceptional current speeds in the channel, regularly exceeding 3.5 m s⁻¹ [45,57,58]. This means that birds drifting passively with the current, a behavior known as tidal drift, could realistically reach groundspeeds normally associated with flight [59,60,61]. An example R script for preparation, processing and analysis (steps 2–4 Figure 2) along with a dataset for one individual are provided in Appendix A and Supplementary Materials.

2.3.2. Utilization Distributions

The extraction of utilization distributions to determine the probability of occurrence of animals is a common method in ecology [62,63]. The outputs are intensity of area use polygons that can be readily interpreted and are highly comparable between sites and species [37]. The 50% and 95% isopleths (i.e., contours) typically extracted from utilization distributions represent the ‘core’ and ‘active’ areas of use, respectively [64]. High-resolution GPS tracks such as those used in this study (positions every 120 s) are autocorrelated in time and space [65,66]. Utilization distributions were therefore extracted using a biased random bridge approach, which accounts for autocorrelation by including a non-random “drift” component to the estimation of the probability density function, implemented in adehabitatHR [52,67]. The diffusion coefficient D (m² s⁻¹, the aggregate of distributions that specify the random-walk model predicting the path of an individual, [68]) was determined using the maximum likelihood function “BRB.likD()” in the adehabitatHR package, described in [69]. The minimum smoothing parameter h_min (m, the minimum uncertainty surrounding an animal’s position) was set to 200 m following the distribution of distances between consecutive fixes [70]. The utilization distributions were computed on a 100 m grid, the size of which was determined by the range of distances between fixes as well as the scales of the tidal draft plan option (>10 km) and tidal lease sites (<5 km). As the purpose of this exercise was to assess potential for overlap with turbine structures at the sea surface (e.g., floating) and underwater (e.g., blades), the ratio of dive locations within the area of use polygons as compared to total number of dives observed was sought to be maximized. The 95% utilization distribution isopleths were therefore extracted for each trip per individual as these retained the majority of dive locations, as determined by the TDRs (Appendix B). Prior to assessment of overlap with tidal sites, areas overlapping land were removed using the sf package in R [71].

2.3.3. Variance and Representativeness

As individuals performing multiple trips may exhibit site fidelity, which in turn may bias results [72,73], a variance test, adapted from [37], was applied. This test determined whether the proximity between areas of use within one individual was significantly different from that of all individuals. First, the proximity (i.e., Hausdorff distances) between 95% occurrence distribution polygons of each trip by individuals that performed more than one trip was calculated [74]. These values were then compared to the reference distribution of distances of all individuals. The reference distribution was attained by randomly selecting an equal number of trips by each individual and calculating proximity between 95% occurrence distributions between individuals. A Mann–Whitney U test was used to compare within individual distances between trips with the reference distribution. The mean p value of 100 iterations of the variance test was calculated. If p > 0.25 then site fidelity was not exhibited and trips from an individual could be used in further analysis, otherwise a random sample trip would need to be taken.

As only a portion of the population is usually tracked, results obtained from tracking studies may not be indicative of population-wide behaviors [75,76,77]. While small tracking sample sizes may be unavoidable due to ethical or logistical considerations [78,79], assessing representativeness can be used as a measure of uncertainty. A test for representativeness was therefore performed using the SDLfilter package [80]. Following methods outlined in [80], the overlap probabilities of the Muckle Skerry utilization distributions, grouped by individual (n = 5) and trip (n = 23), were quantified, as per the recommendations of [81]. For each group, the function “boot_overlap()” was applied, with 10,000 iterations and method “PHR”(the probability distribution) specified. The overlap probability calculates the probability of a sample utilization distribution being within the utilization distribution of another sample, thereby quantifying the contribution of each new UD to the existing samples. The function calculates mean overlap probability for each sample size following built-in boot-strapping and estimates a horizontal asymptote and curve approaching it using a rational function. The value at which 95% of the horizontal asymptote is crossed indicates the minimum sample size needed to represent the overall distribution (i.e., representativeness).

2.4. Overlap with Tidal Development and Lease Sites

The two-dimensional percentage overlap of the 95% utilization distribution isopleths with the tidal draft plan option and tidal lease sites (e.g., percent of site covered by the isopleth) within the Pentland Firth were quantified per individual and for all the shags using the sf package in R [71].

3. Results

Summary metrics of the at-sea shag data used in the utilization distribution analysis are presented in Table 1. Variance between individuals was found not to be significant (p = 0.26), meaning that tracks from all individuals could be used in subsequent analysis. Greater than 95% overlap probability was achieved at both individual and trip levels, indicating adequate representativeness (Appendix C). The entire shag foraging occurrence was contained within the tidal draft plan option (99.4%), and this accounted for 4.6% of the option area. There was minimal overlap between shag foraging occurrence and the tidal lease sites (0.14%), of which 0.5% overlapped with the foraging occurrence of individual EUSH_616 (

Table 2. Size of the 95% isopleths of foraging distributions per individual and combined for all shags, and percentage overlaps with the tidal draft plan option site (Pentland Firth) and tidal lease sites within (Inner Sound and Ness of Duncansby).

ID	Area (km2)	Tidal PO Overlap		Tidal Lease Site Overlap
		% of Shag UD	% of Tidal PO	% of Shag UD	% of Tidal Lease Site
EUSH616	7.7	99.9	1.7	0.4	0.5
ORK0445	10.3	99.0	2.3	0	0
ORK0451	9.9	99.8	2.2	0	0
ORK0694	0.6	100.0	0.13	0	0
ORK0795	3.5	100.0	0.78	0	0
ORK0797	2.6	100.0	0.59	0	0
Total	20.5	99.4	4.6	0.14	0.5

, Figure 3).

4. Discussion

The potential adverse effects of marine renewable energy devices on marine wildlife, for example the risk of underwater collision with tidal turbines, remain uncertain [5]. In addition, the need to assess the magnitude of any effects and adopt a precautionary approach is likely to delay or even prevent developments where empirical measurements of overlap between animals and sites are lacking. This could make it more difficult to reach net-zero emissions [2,3,83]. The potential for animal-borne devices to reveal areas of use and therefore potential for overlap with marine renewables is great [26,84,85,86,87]. However, there are ethical and logistical concerns associated with tagging wild animals, as well as time and resource constraints to consider when analyzing the complex data retrieved from the devices [17,41,79]. Additionally, while such data have been used in marine conservation and management previously (e.g., [26,85]), this is often at a scale (i.e., ocean-basin) less immediately relevant to developers and stakeholders in marine renewables. In order to facilitate its use within the assessment of risk to animals from marine renewable energy, we developed an analysis workflow (Figure 2) and applied it to archival biologging and telemetry data for six European shags within a major tidal development area.

Within the field of animal-borne biologging and telemetry, there is a broad range of data preparation, processing, and analysis techniques, making comparability and reproducibility challenging [87,88,89,90,91]. For instance, the sheer number of R packages to analyze animal movement alone may well be overwhelming to new practitioners [92]. There is therefore a need for clear and reproducible workflows to facilitate deriving area of use polygons for incorporation in marine environmental impact assessments or future siting of developments, similar to that developed for the identification of important seabird conservation sites from GPS-only data (e.g., [37,93,94,95]).

The framework presented here expands on Lascelles et al. [37] in that the stressor of interest (e.g., collision) influences choices made throughout (steps 2–4 in Figure 2), relevant biologging data are incorporated (steps 2–4 in Figure 2) and there is greater flexibility in the choice of utilization distribution, isopleths extracted, and calculation of representativeness (step 4 in Figure 2). Furthermore, the preparation and processing steps required are made explicit, similar to the approach in appendix 2 of Handley et al. [27]. It is intended as a malleable guide, allowing the user to modify specific aspects (dotted boxes in Figure 2) depending on the stressor, species, tags, sampling regime, and familiarity with analytical methods. It is worth noting that for each step in the workflow there are a myriad of packages available in R (reviewed in [92]) and that the worked example given in this study (Appendix A) identifies a few relevant ones. In the interest of accessibility and reproducibility, R was the software of choice (R is a major coding language in movement ecology), yet other open-source programs such as Python or QGIS may be alternatives.

Several of the components of the analysis workflow merit further discussion. The amount of processing in step two (data preparation, Figure 2) will vary depending on the type of animal-borne device deployed, sampling regime and time period sampled. Visual inspection for outliers and removal of duplicates, and, where available, filtration based on metrics of positional quality often recorded in conjunction with GPS fixes, such as dilution of precision or number of satellites, should be performed. Seabird tracking data are often collected during the breeding season, when the birds are central-place foragers and therefore necessitating analysis at a trip level, as performed here [53]. This processing step may not be appropriate for data collected outside of this season nor for other marine animals. Processing in step three (data processing, Figure 2) will depend on the species and behavior of interest which in turn depends on the stressor applied by the renewable energy technology. In this instance, only foraging positions (on the sea surface and below) were retained due to shags being pursuit-divers and the energy extracted being tidal [96]. For example, in the case of plunge-divers such as Northern Gannet (Morus bassanus) or overlap with offshore wind, positions identified as flight might be retained [97]. There are also many methods to identify animal behaviors including first passage time, speed/tortuosity thresholds, state-space models and machine learning [48,56,98,99,100]. As dive locations were known from the TDRs, a speed and tortuosity threshold was applied to distinguish flight from foraging, although any of the other methods mentioned previously could be substituted in the workflow (step 3, Figure 2).

In order to extract information on space use from the processed points, utilization distributions, commonly used in ecology to identify areas used by animals based on positional data, were derived [62,64] (step 4, Figure 2). Utilization distributions are used to estimate where an animal is located over a season (home range) or over a short tag deployment period (occurrence distribution) [65]. Due to the short sampling period (on the scale of days, see Table 1) per individual in the shag example, the utilization distribution calculated represented an occurrence distribution [47,65]. As tracking data are serially autocorrelated (points near each other in time tend to be near each other in space as well) at high sampling intervals, a biased random bridge estimator was used although other estimators such as Continuous-time correlated RAndom Walk (CRAWL) or timeseries Kriging may be appropriate [67,70,101,102]. Where data are taken to be representative of the individual’s use of space over a longer timespan (e.g., an entire breeding season), a home range estimator such as kernel density may be applied instead, which would likely yield larger isopleths [62]. Larger isopleths are more precautionary and may therefore be more appropriate in impact assessments; the use of a longer timeseries tracking dataset and autocorrelated kernel density estimator (AKDE), appropriate for high-resolution tracking data, would achieve this [65,103]. The workflow allows for flexibility in choice of distribution, and estimator, depending on data, level of caution, and preference of the practitioner (see step 4, Figure 2).

The isopleths commonly extracted from utilization distributions are typically taken to represent the ‘core’ (50%) and ‘active’ (95%) areas [64,81]. Whether core or active areas should be used in further analysis will depend a great deal on the level of caution being ascribed, which in turn depends on the stressor and species. In the example given, where shag collision with turbine structures either at the surface (e.g., floating mooring) and especially underwater (e.g., rotating blades) is highly undesirable, the 95% isopleths were extracted (step 4, Figure 2). Comparing the ratios of dives retained within the 50 and 95% isopleths revealed that a majority of known dive locations (>94%) were retained within the 95% isopleth, while anywhere from 18 to 86% were retained in the 50% isopleth (Appendix B). In the more common case where locations of dives are not known (i.e., GPS-only tracking data), and given similar behavioral analyses, this result can aide isopleth level selection.

Within animal-borne tracking datasets, individuals often contribute data unequally. For example, over the same time period, an individual may perform more trips while another from the same tagging campaign remains relatively stationary. This may bias results should individuals exhibit site fidelity [72,73]. A test to determine whether all trips by a central-place foraging individual can be used in subsequent analyses or whether a random sample needs to be taken should therefore be applied [37] (step 4, Figure 2). The extremely conservative significance threshold (p = 0.25) recommended by Lascelles et al. [37] will depend on the species’ overall tendency to site fidelity, although should be set no lower than 0.05. As animal-borne tracking data are by definition obtained for individuals, it is not necessarily representative of the wider population and therefore a test for representativeness is also recommended [75,81,94] (step 4, Figure 2, Appendix C). Small sample sizes may be unavoidable, yet in this way the level of representativeness can serve as a proxy for uncertainty. This can inform risk level applied to the resulting areas of use as well as determine the need and extent of future data collection efforts.

The area of use polygons resulting from the workflow have applications beyond the assessment of the extent of spatial overlap animals have with marine renewable energy developments. The polygons can be used by stakeholders to determine connectivity between a colony and a renewable energy development area as well as flag exposure to collision risk at specific sites. This is useful as it identifies populations that a development project may impact and, if a Special Protected Area (SPA) population, indicate need for an Appropriate Assessment under the EU Birds or Habitats Directives [104,105]. The representativeness of the polygons could also help to determine the need for additional data to increase sample size and/or temporal coverage (e.g., seasons, years). Should the polygons be deemed sufficiently robust, relevant complementary data from devices (e.g., depth) could be used to populate collision risk models to estimate collision mortality [106,107,108,109,110].

In the worked example in this study, the difference in overlap between shags and the draft plan option and tidal lease sites in the Pentland Firth highlights the importance of scale [12,13,111] (Table 2, Figure 3). At the scale of >10 km, the overlap between seabird area use and development area was 99.4%, whereas at the much smaller scale at which tidal lease sites occur in the UK (<5 km), there was negligible overlap. It is important to note here that while the data were assessed to be adequately representative of Muckle Skerry shag occurrence, this is explicitly for a few days within the breeding season based on data collected prior to installation of turbines within the site (Appendix C). Shags are also fairly local foragers (maximum range ca. 17 km, [82]), and the degree of overlap will almost certainly depend on proximity of sites to the breeding colonies, (Figure 3B). Shags nesting along the mainland Scottish coast just south of the Inner Sound or Ness of Duncansby lease sites should therefore be prioritized in data collection efforts (Figure 3B,C). Incidentally, the directional shape of the distribution suggests that the shags favour areas conducive to foraging (rather than uniformly foraging within a 17 km radius, Figure 3A), that may readily be explained by environmental factors such as bathymetry, hydrodynamics, or prey [13,112,113]. It is also worth noting that the same workflow performed on diving species that forage farther afield (e.g., razorbills Alca torda, harbor seals Phoca vitulina) will likely yield larger polygons. While density of birds within the distribution is not explored further here (some indication is given in Figure 3C), this becomes increasingly relevant when considering non-uniform placement of arrays of energy extracting devices within the development area.