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

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.


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].
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 −1 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-windwave-and-tidal-2013, accessed on 1 July 2020 [14]) and Crown Estate Scotland (https: //www.crownestatescotland.com/, accessed on 1 February 2020).

Bird-Borne Biologging and Telemetry Dataset
Data from European shags simultaneously carrying GPS tags and TDRs tained through the Environmental Research Institute and RSPB FAME/STAR (www.fameproject.eu, accessed on 1 October 2018). The dataset consisted of tagged within the Pentland Firth during the breeding season (May-June; Five o Skerry (2012 n = 2; 2013 n = 1; 2014 n = 2) and one on Stroma (2012 n = 1). (Figure 1 Appendix A). The combined weight of devices was always <3% of shag body we further details on bird handling, attachment method and devices can be found in Table 1. Summary statistics of the six shags included in the analysis.

Preparation and Processing
Data preparation, processing and analysis were carried out in R version 3. Studio version 1.2.1335 [51]. The steps for preparation and processing, adapted are outlined in Figure 2 and followed methods outlined in [48]. The simultan and TDR datasets from each individual were merged to the nearest GPS poin resulting in a dataset with dive locations identified. The dataset was then red (linearly interpolated) to 120 s intervals between fixes using the adehabitatLT pac

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.  Table 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].

. 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].

Figure 2.
Analysis workflow to go from raw bird-borne biologging and telemetry point data to area of use po assess overlap with marine renewable sites. After obtaining animal-borne biologging and telemetry data in the the development area (step 1), preparation (step 2) consists mainly of quality control before identification of relevant to the stressor being assessed (step 3), and extraction of utilization distributions (step 4). The workflow as a guide and specifics within each step will depend on characteristics of the data, species, and renewable en nology.
During the breeding season, shags are a central place foragers and are strained to a maximum distance from the nest, making trips from the ne foraging sites in order to defend, incubate and feed their offspring [53]. For sification is therefore an inherent part of breeding seabird track analysis. were defined as a minimum 30 min round trip from and to a 300 m buffer ar These thresholds were chosen to eliminate resting and washing behaviors ble [54,55]. As the stressor of interest was collision with tidal energy dev deployed at, in the case of floating tidal, or below the sea surface, GPS fix flight as well as those on land were removed, leaving only at-sea fixes like foraging (dives or sitting on the sea surface). Flight fixes were identified ba and tortuosity threshold of groundspeeds over 5 m s −1 and turning angles o [55,56]. The speed threshold was chosen to account for the exceptional cur the channel, regularly exceeding 3.5 m s −1 [45,57,58]. This means that bird sively with the current, a behavior known as tidal drift, could realistically speeds normally associated with flight [59][60][61]. An example R script for pre cessing and analysis (steps 2-4 Figure 2) along with a dataset for one indiv vided in Appendix A and Supplementary Materials.

Utilization Distributions
The extraction of utilization distributions to determine the probability of animals is a common method in ecology [62,63]. The outputs are intens polygons that can be readily interpreted and are highly comparable betw species [37]. The 50% and 95% isopleths (i.e., contours) typically extracted fr distributions represent the 'core' and 'active' areas of use, respectively [64 tion GPS tracks such as those used in this study (positions every 120 s) are . Analysis workflow to go from raw bird-borne biologging and telemetry point data to area of use polygons and assess overlap with marine renewable sites. After obtaining animal-borne biologging and telemetry data in the vicinity of the development area (step 1), preparation (step 2) consists mainly of quality control before identification of behaviors relevant to the stressor being assessed (step 3), and extraction of utilization distributions (step 4). The workflow is intended as a guide and specifics within each step will depend on characteristics of the data, species, and renewable energy technology.
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 −1 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 −1 [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.

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 2 s −1 , 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].

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 bootstrapping 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).

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].

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, Figure 3). = 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, Figure 3).

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. Funding: This work was funded by the Bryden Centre project, supported by the European Union's INTERREG VA Programme, managed by the Special EU Programmes Body (SEUPB). The APC was funded by PRIMaRE. FAME and STAR Funding in Orkney was provided by the European Regional Development Fund through its Atlantic Area Programme, Marine Scotland, Scottish Natural Heritage, the Joint Nature Conservation Committee and the RSPB.

Institutional Review Board Statement:
This study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Research Ethics Committee of University of the Highlands and Islands (OL-ETH SHE-1230, approved on 21 February 2019). All RSPB logger deployments were carried out following the ethical guidelines of the British Trust for Ornithology, under license by Scottish Natural Heritage (for more details see: [113]).

Informed Consent Statement: Not applicable.
Data Availability Statement: Data from five shags tagged on Muckle Skerry available on request to RSPB (https://www.rspb.org.uk/our-work/conservation/projects/tracking-seabirds-to-informconservation-of-the-marine-environment/, accessed on 31 December 2020). Data from one shag tagged on Stroma is provided within Supplementary Material. and/or drafts. We gratefully acknowledge the many fieldworkers who assisted with data collection (Derren Fox, Neil James, Juliet Lamb, Kirsty Lees, Tegan Newman, Yvan Satge, Chris Taylor, Jenni Taylor, and others). Thanks to three anonymous reviewers for comments that improved this paper.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Worked Example of Data Preparation, Processing and Analysis for One Shag
The following script details steps 2-4 (excluding variance and representativeness analysis) in the workflow for the use of simultaneous animal-borne biologging and telemetry data to arrive at utilization distribution areas for use in assessing overlap with marine renewables (Figure 2). The following worked example uses data from a shag tagged on Stroma in 2012 'named' EUSH_616, raw data files provided in Supplementary Materials. and/or drafts. We gratefully acknowledge the many fieldworkers who assisted with data collection (Derren Fox, Neil James, Juliet Lamb, Kirsty Lees, Tegan Newman, Yvan Satge, Chris Taylor, Jenni Taylor, and others). Thanks to three anonymous reviewers for comments that improved this paper.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Worked Example of Data Preparation, Processing and Analysis for One Shag
The following script details steps 2-4 (excluding variance and representativeness analysis) in the workflow for the use of simultaneous animal-borne biologging and telemetry data to arrive at utilization distribution areas for use in assessing overlap with marine renewables (Figure 2). The following worked example uses data from a shag tagged on Stroma in 2012 'named' EUSH_616, raw data files provided in Supplementary Materials.

Read in TDR and GPS Data
Before beginning, the TDR dataset was 'clipped' to the same deployment length as the GPS, and pressure readings (in dBar) were converted to depth in the water column (m) using a conversion factor of 1 Bar = 10.1974 m.

Create TDR Object
Create a TDR object for working in divemove package [114,115]. Check the metadata and plot the dives. and/or drafts. We gratefully acknowledge the many fieldworkers who assisted with data collection (Derren Fox, Neil James, Juliet Lamb, Kirsty Lees, Tegan Newman, Yvan Satge, Chris Taylor, Jenni Taylor, and others). Thanks to three anonymous reviewers for comments that improved this paper.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Worked Example of Data Preparation, Processing and Analysis for One Shag
The following script details steps 2-4 (excluding variance and representativeness analysis) in the workflow for the use of simultaneous animal-borne biologging and telemetry data to arrive at utilization distribution areas for use in assessing overlap with marine renewables (Figure 2). The following worked example uses data from a shag tagged on Stroma in 2012 'named' EUSH_616, raw data files provided in Supplementary Materials.

Read in TDR and GPS Data
Before beginning, the TDR dataset was 'clipped' to the same deployment length as the GPS, and pressure readings (in dBar) were converted to depth in the water column (m) using a conversion factor of 1 Bar = 10.1974 m.

Choose Smoothing Parameters
After visual inspection of data and with reference to package documentation, select smoothing parameters. These will allow for the surface to be set and for drift to be corrected. k = c(x, y), where x is the window width (in seconds) of the first smoothing step (median), i.e., the small one, (in this case 3 s given 1 s intervals), and y is the window width for the large filter (in this case 120 s). p = c(a, b), where a is the quantile for the first step and b is the quantile for the second, if there is little noise after first filter then a smaller quantile (i.e., 0.01-0.05) can be used for the second step.

Choose Smoothing Parameters
After visual inspection of data and with reference to package documentation, select smoothing parameters. These will allow for the surface to be set and for drift to be corrected. k = c(x, y), where x is the window width (in seconds) of the first smoothing step (median), i.e., the small one, (in this case 3 s given 1 s intervals), and y is the window width for the large filter (in this case 120 s). p = c(a, b), where a is the quantile for the first step and b is the quantile for the second, if there is little noise after first filter then a smaller quantile (i.e., 0.01-0.05) can be used for the second step.

Calibrate Data
Select k and p values from previous step and extract dives with calibrated data. In this case, dive threshold was set to 3 m for consistency with method in [48]. After visual inspection of data and with reference to package documentation, select smoothing parameters. These will allow for the surface to be set and for drift to be corrected. k = c(x, y), where x is the window width (in seconds) of the first smoothing step (median), i.e., the small one, (in this case 3 s given 1 s intervals), and y is the window width for the large filter (in this case 120 s). p = c(a, b), where a is the quantile for the first step and b is the quantile for the second, if there is little noise after first filter then a smaller quantile (i.e., 0.01-0.05) can be used for the second step.
old.par<-par(no.readonly = TRUE) plotZOC(shag, d.filter, ylim =c(-1.029937, 5)) Visualize smoothing parameters. First row is original data, second is smoothing, and 3 is the what the data looks like with the smoothing applied.

Choose Smoothing Parameters
After visual inspection of data and with reference to package documentation, select smoothing parameters. These will allow for the surface to be set and for drift to be corrected. k = c(x, y), where x is the window width (in seconds) of the first smoothing step (median), i.e., the small one, (in this case 3 s given 1 s intervals), and y is the window width for the large filter (in this case 120 s). p = c(a, b), where a is the quantile for the first step and b is the quantile for the second, if there is little noise after first filter then a smaller quantile (i.e., 0.01-0.05) can be used for the second step.

Calibrate Data
Select k and p values from previous step and extract dives with calibrated data. In this case, dive threshold was set to 3 m for consistency with method in [48].

Calibrate Data
Select k and p values from previous step and extract dives with calibrated data. In this case, dive threshold was set to 3 m for consistency with method in [48].

Calibrate Data
Select k and p values from previous step and extract dives with calibrated data. In this case, dive threshold was set to 3 m for consistency with method in [48].

# Function to get the index specifying closest or after
Merge the gps and tdr datasets. #attr(tdr_final$merge_date, "tzone") #"" tz(gps$merge_date)<-"UTC" #f not identical or timezone not assigned, assign tz(tdr_final$merge_date)<-"UTC" #f not identical or timezone not assigned, assign tdr_final$merge_date<-round_date(tdr_final$merge_date, "minute") #round to nearest minute to match gps df Merge the gps and tdr datasets. Now that the raw biologging data has been processed and merged with the GPS data, resulting in a dataset with known dive locations, the track should be visualized ( Figure  2). Base maps should be high enough resolution for subsequent analyses (e.g., removal of land), so in this case were downloaded from Digimap: https://digimap.edina.ac.uk/os (accessed on 1 February 2020). There are no obvious outliers or positions extremely far from the colony.

Tidy
Tidy up resulting dataframe of 587 observations and 79 dives.
#attr(tdr_final$merge_date, "tzone") #"" tz(gps$merge_date)<-"UTC" #f not identical or timezone not assigned, assign tz(tdr_final$merge_date)<-"UTC" #f not identical or timezone not assigned, assign tdr_final$merge_date<-round_date(tdr_final$merge_date, "minute") #round to nearest minute to match gps df Merge the gps and tdr datasets. Now that the raw biologging data has been processed and merged with the GPS data, resulting in a dataset with known dive locations, the track should be visualized ( Figure  2). Base maps should be high enough resolution for subsequent analyses (e.g., removal of land), so in this case were downloaded from Digimap: https://digimap.edina.ac.uk/os (accessed on 1 February 2020). There are no obvious outliers or positions extremely far from the colony.

Appendix A.1.2. Visualize
Now that the raw biologging data has been processed and merged with the GPS data, resulting in a dataset with known dive locations, the track should be visualized ( Figure 2). Base maps should be high enough resolution for subsequent analyses (e.g., removal of land), so in this case were downloaded from Digimap: https://digimap.edina.ac.uk/os (accessed on 1 February 2020). There are no obvious outliers or positions extremely far from the colony. When these variables are available, data can be filtered to only include points with low dilution of precision and/or high number of satellites. This is not the case in this example, which instead will require paying close attention to groundspeeds calculated further down in the analysis and removing any fixes that reach above maximum flight speeds of shags (>30 m/s [117]).

Appendix A.2. Step 3: Data Processing
Appendix A.

Rediscretize to Consistent Sampling Interval
The package adehabitatLT [52] can be used to regularize the sampling interval, in this case to 120 s between each fix. Converting the dataframe into an 'ltraj' object required for the rediscretization function to work will calculate a number of relevant variables, including distance covered between succesive fixes ('dist'). Should a readily interpretable unit (e.g., metres) be desired, the coordinate reference system (CRS) supplied needs to be in that unit. In this case, the CRS is WGS84 (standard for GPS data), so Longitude and Latitude are in decidegrees. Prior to rediscretization, this can be converted to a CRS with units in metres, British National Grid in this example, although Universal Transverse Mercator would also work. The 'st_transform' and 'map' functions contained in the sf and tidyverse packages can be used for this conversion [71,118].

Remove Low Quality Positions
When these variables are available, data can be filtered to only include points with low dilution of precision and/or high number of satellites. This is not the case in this example, which instead will require paying close attention to groundspeeds calculated further down in the analysis and removing any fixes that reach above maximum flight speeds of shags (>30 m/s [117]). The package adehabitatLT [52] can be used to regularize the sampling interval, in this case to 120 s between each fix. Converting the dataframe into an 'ltraj' object required for the rediscretization function to work will calculate a number of relevant variables, including distance covered between succesive fixes ('dist'). Should a readily interpretable unit (e.g., metres) be desired, the coordinate reference system (CRS) supplied needs to be in that unit. In this case, the CRS is WGS84 (standard for GPS data), so Longitude and Latitude are in decidegrees. Prior to rediscretization, this can be converted to a CRS with units in metres, British National Grid in this example, although Universal Transverse Mercator would also work. The 'st_transform' and 'map' functions contained in the sf and tidyverse packages can be used for this conversion [71,118]. When these variables are available, data can be filtered to only include points with low dilution of precision and/or high number of satellites. This is not the case in this example, which instead will require paying close attention to groundspeeds calculated further down in the analysis and removing any fixes that reach above maximum flight speeds of shags (>30 m/s [117]).

Appendix A.2. Step 3: Data Processing
Appendix A.

Rediscretize to Consistent Sampling Interval
The package adehabitatLT [52] can be used to regularize the sampling interval, in this case to 120 s between each fix. Converting the dataframe into an 'ltraj' object required for the rediscretization function to work will calculate a number of relevant variables, including distance covered between succesive fixes ('dist'). Should a readily interpretable unit (e.g., metres) be desired, the coordinate reference system (CRS) supplied needs to be in that unit. In this case, the CRS is WGS84 (standard for GPS data), so Longitude and Latitude are in decidegrees. Prior to rediscretization, this can be converted to a CRS with units in metres, British National Grid in this example, although Universal Transverse Mercator would also work. The 'st_transform' and 'map' functions contained in the sf and tidyverse packages can be used for this conversion [71,118].

Appendix A.2.2. Set Colony/Nest Buffer and Assign Trips
During breeding, shags are a central place foragers and perform regular foraging trips from and to the nest [53]. In order to identify trips, fixes can be assigned as being: 0: not on trip; 1: last point of trip; 3: middle points of trip; 2: first point of trip; where the points are relative to a biologically relevant buffer distance set around the nest. In this case, 300 m was chosen in order to eliminate resting and washing behaviors as far as possible [54,55,119]. First, distance to the colony needs to be calculated.

Set Colony/Nest Buffer and Assign Trips
During breeding, shags are a central place foragers and perform regular foraging trips from and to the nest [53]. In order to identify trips, fixes can be assigned as being: 0: not on trip; 1: last point of trip; 3: middle points of trip; 2: first point of trip; where the points are relative to a biologically relevant buffer distance set around the nest. In this case, 300 m was chosen in order to eliminate resting and washing behaviors as far as possible [54,55,119]. First, distance to the colony needs to be calculated.  (activity), "D", as.character(activity))) The rediscretized dataset now contains 698 observations.

Appendix A.2.2. Set Colony/Nest Buffer and Assign Trips
During breeding, shags are a central place foragers and perform regular foraging trips from and to the nest [53]. In order to identify trips, fixes can be assigned as being: 0: not on trip; 1: last point of trip; 3: middle points of trip; 2: first point of trip; where the points are relative to a biologically relevant buffer distance set around the nest. In this case, 300 m was chosen in order to eliminate resting and washing behaviors as far as possible [54,55,119]. First, distance to the colony needs to be calculated. #add nest location to df as_sf2$nest_geom<-st_as_sf(interp[,c("nest_x","nest_y")], coords = c("nest_x","nest_y"),crs=27700) #compute colony distance variable as_sf2$col_dist <-st_distance(as_sf2$geometry, as_sf2$nest_geom, by_element = TRUE) #default Euclidean Now fixes can be assigned numbers based on where they fall in relation to the buffer, and from here assigned to a trip.
The following function (trip_assign) can do this, provided there is a (numeric) 'col_dist' and 'id' variable in the dataframe.
The following function (trip_assign) can do this, provided there is a (numeric) 'col_dist' and 'id' variable in the dataframe.

(i in 2:nrow(df)){ #go through df starting from the second row if(df$id[i]==df$id[i-1]){ #check if animal id is the same as the animal id of the row before if(dist_over[i]==1){ #check if the animal is on a trip if(dist_over[i]==dist_over[i-1]){ #check if it was on a trip the previous row trip_id[i]<-current_trip #assign the number of the current trip in the list } else{ #if animal is on a new trip current_trip<-current_trip+1 #assign trip +1 trip_id[i]<-current_trip } }else{ #if animal is not on a trip trip_id[i]<-0 } }else{ #if its a new animal, start from 1 or 0 again depending on if the animal is on a trip or not current_trip<-dist_over[i] trip_id[i]<-current_trip } } return(trip_id) }
#make sure df is ordered! as_sf2 <-as_sf2[order(as_sf2$id, as_sf2$date),] #df must be ordered by id and consecutive fixes as_sf2$trip_id<-trip_assign(as_sf2, 300) #assign trip based on 300m buffer using trip_assign() provided #anything with trip id '0' is within buffer so remove EUSH_616_trips<-as_sf2%>%filter(trip_id != 0) #in some cases it may be of interest to see how much diving occurs within the buffer zone #EUSH_616_nontrip<-as_sf2%>%filter(trip_id == 0) #(sum(nrow(EUSH_616_nontrip %>% filter(activity == "W")))/nrow(EUSH_616_nontrip))*100 #0.27% Applying a duration (e.g., >30 min duration) and minimum maximum distance (e.g., 500 m) to nest filtration on trips helps ensure likely foraging behavior as well [54,55]. This individual performed one trip with 324 fixes, with a duration of 10.7 h and a maximum distance of 16.1 km from the nest. In this case, no further filtration on the trip is required. Applying a duration (e.g., >30 min duration) and minimum maximum distance (e.g., 500 m) to nest filtration on trips helps ensure likely foraging behavior as well [54,55]. This individual performed one trip with 324 fixes, with a duration of 10.7 h and a maximum distance of 16 The dataset now needs to be filtered to exclude points that are likely flight and on land. This is because for the purpose of further analysis we are only interested in observations either on the surface of the sea or below (where components of tidal turbines are).
One straightforward method to distinguish between flight and rest/foraging is to apply both speed and tortuosity thresholds [55,56,[120][121][122]. This works on the principle that flight is associated with high speed and low tortuosity. In other words, commuting/transit tends to have a high degree of directionality, that is, tends to occur in a straight line. Foraging is associated with lower speed and high tortuosity, and resting with low speed and low tortuosity [55,56,123].

Appendix A.3. Speed
In order to identify flight fixes, groundspeed needs to be calculated. As difference in distance and time variables have previously been calculated when converting to an ltraj object, this is simply a matter of dividing 'dist' by 'dt'.

Appendix A.2.3. Identify Behaviors of Interest
The dataset now needs to be filtered to exclude points that are likely flight and on land. This is because for the purpose of further analysis we are only interested in observations either on the surface of the sea or below (where components of tidal turbines are).
One straightforward method to distinguish between flight and rest/foraging is to apply both speed and tortuosity thresholds [55,56,[120][121][122]. This works on the principle that flight is associated with high speed and low tortuosity. In other words, commuting/transit tends to have a high degree of directionality, that is, tends to occur in a straight line. Foraging is associated with lower speed and high tortuosity, and resting with low speed and low tortuosity [55,56,123].

Appendix A.3. Speed
In order to identify flight fixes, groundspeed needs to be calculated. As difference in distance and time variables have previously been calculated when converting to an ltraj object, this is simply a matter of dividing 'dist' by 'dt'. Examining a histogram of these ground speeds reveals a range that agrees with those reported for shags in the literature [96,117,124].
The minimum groundspeed of a flying shag is 1 m/s [117], shown in the above plot with a dashed line. The blue shaded area (3.5-5 m/s) is the range of peak current speeds in the Pentland Firth [45,57,58]. Between 1 and 5 m/s, there is therefore a range of speeds Figure A5. Speed histogram.
Examining a histogram of these ground speeds reveals a range that agrees with those reported for shags in the literature [96,117,124].
The minimum groundspeed of a flying shag is 1 m/s [117], shown in the above plot with a dashed line. The blue shaded area (3.5-5 m/s) is the range of peak current speeds in the Pentland Firth [45,57,58]. Between 1 and 5 m/s, there is therefore a range of speeds that normally could be classified as flight but in this case could be tidal drift (bird sitting on the sea surface and drifting with the tide, [59][60][61]). The cautious approach would therefore be to set a speed threshold of >5 m/s.

Appendix A.4. Tortuosity
The 'rel.angle' (relative angle) variable computed by adehabitatLT used when rediscretizing the tracks to 120 s measures the change of direction between the step built by fix − 1 and the current fix and the step between current fix and fix + 1, also known as the turning angle [52].
In order for the turning angle to be readily interpretable and comparable with other values in the literature, we can calculate a tortuosity index that goes from 0 (tortuous) to 1 (straight). Examining a histogram of these ground speeds reveals a range that agrees with those reported for shags in the literature [96,117,124].
The minimum groundspeed of a flying shag is 1 m/s [117], shown in the above plot with a dashed line. The blue shaded area (3.5-5 m/s) is the range of peak current speeds in the Pentland Firth [45,57,58]. Between 1 and 5 m/s, there is therefore a range of speeds that normally could be classified as flight but in this case could be tidal drift (bird sitting on the sea surface and drifting with the tide, [59][60][61]). The cautious approach would therefore be to set a speed threshold of >5 m/s.

Appendix A.4. Tortuosity
The 'rel.angle' (relative angle) variable computed by adehabitatLT used when rediscretizing the tracks to 120 s measures the change of direction between the step built by fix − 1 and the current fix and the step between current fix and fix + 1, also known as the turning angle [52].
In order for the turning angle to be readily interpretable and comparable with other values in the literature, we can calculate a tortuosity index that goes from 0 (tortuous) to 1 (straight). The red dashed line shows a turning angle equivalent to 45 degrees, previously set in shag behavioral classification analysis [55], above which the fix is assessed to be straight.

Appendix A.5. Identify Flight Using Speed, Tortuosity and Dive Classification
Now that groundspeed and a tortuosity index have been calculated, these can be combined with the dive locations identified via TDR to determine behavioral states and identify flight. The red dashed line shows a turning angle equivalent to 45 degrees, previously set in shag behavioral classification analysis [55], above which the fix is assessed to be straight.  The behaviors that can be discriminated in the above plot are resting (bottom right quadrant), foraging (bottom and top left quadrants), and transit/flight (top right quadrant). Where fixes are marked as dives these take precedence over speed and tortuosity classifications, as these are identified empirically. For example, 1 dive fix occurs at transit speeds and this is likely due to the fast tidal currents in the area.

Retain Only Behaviors of Interest
The dataset can now be filtered by dive classification, speed, and tortuosity to remove flight, leaving 297 observations.

noland<-st_difference(EUSH_616_nonflight, scot)
The resulting dataset contains 236 observations. The behaviors that can be discriminated in the above plot are resting (bottom right quadrant), foraging (bottom and top left quadrants), and transit/flight (top right quadrant). Where fixes are marked as dives these take precedence over speed and tortuosity classifications, as these are identified empirically. For example, 1 dive fix occurs at transit speeds and this is likely due to the fast tidal currents in the area.

Retain Only Behaviors of Interest
The dataset can now be filtered by dive classification, speed, and tortuosity to remove flight, leaving 297 observations. The behaviors that can be discriminated in the above plot are resting (bottom right quadrant), foraging (bottom and top left quadrants), and transit/flight (top right quadrant). Where fixes are marked as dives these take precedence over speed and tortuosity classifications, as these are identified empirically. For example, 1 dive fix occurs at transit speeds and this is likely due to the fast tidal currents in the area.

Retain Only Behaviors of Interest
The dataset can now be filtered by dive classification, speed, and tortuosity to remove flight, leaving 297 observations. EUSH_616_nonflight<-EUSH_616_trips%>%filter(activity == "W"| speed_ms<5 | t_index <0.7) #297 obs The remaining dataset can then be overlaid on a land base map and sites on land (e.g., roosts) can then be removed using the st_difference function in the sf package, leaving only foraging/resting-at-sea points [71].
The remaining dataset can then be overlaid on a land base map and sites on land (e.g., roosts) can then be removed using the st_difference function in the sf package, leaving only foraging/resting-at-sea points [71]. The behaviors that can be discriminated in the above plot are resting (bottom right quadrant), foraging (bottom and top left quadrants), and transit/flight (top right quadrant). Where fixes are marked as dives these take precedence over speed and tortuosity classifications, as these are identified empirically. For example, 1 dive fix occurs at transit speeds and this is likely due to the fast tidal currents in the area.

Retain Only Behaviors of Interest
The dataset can now be filtered by dive classification, speed, and tortuosity to remove flight, leaving 297 observations. EUSH_616_nonflight<-EUSH_616_trips%>%filter(activity == "W"| speed_ms<5 | t_index <0.7) #297 obs The remaining dataset can then be overlaid on a land base map and sites on land (e.g., roosts) can then be removed using the st_difference function in the sf package, leaving only foraging/resting-at-sea points [71].

. Step 4: Analysis
Now the track has been prepared and processed to only include foraging locations (positions at which the shag is either on the sea surface or below). This is where spatial overlap with either floating mooring or rotating blades of tidal stream turbines is possible. The next step is extract the utilization distribution (UD) from this track, effectively transforming the point data into a volume. As the length of the track is short (<48 h), the UD will be one of occurrence as opposed to a home range [62,64]. As the track is autocorrelated in time and space [65,66], the biased random bridge method implemented in the adehabitatHR package was applied [52]. Where data are taken to be representative of the individual's use of the space over a longer timespan (e.g., an entire breeding season), a home range estimator such as auto-correlated kernel density estimation may be applied instead, implemented in the 'ctmm' package [103,125]. Appendix A.6.1. Utilization Distribution The "BRB()" function in adehabitatHR accounts for autocorrelation by including a non-random "drift" component to the estimation of the probability density function [67]. The diffusion coefficient D (in units of the coordinates, in this case m 2 /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()", described in more detail in [69]. The minimum smoothing parameter hmin (m, the minimum uncertainty surrounding an animal's position) was set to 200 m following the distribution of distances between consecutive fixes [70]. Grid size will depend on study species' movement ecology (range of distances between fixes as well as maximum range), GPS interval and error, scale of environmental variables being compared with, and computing power. The smaller the grid size, the higher the resolution, the larger the file. In this case, grid size (100 m) was selected based on range of distances between fixes, scale of tidal lease site (<5 km), and reasonable file size.
ng<-noland a 616_foraging$x) 616_foraging$y) Now the track has been prepared and processed to only include foraging locations (positions at which the shag is either on the sea surface or below). This is where spatial overlap with either floating mooring or rotating blades of tidal stream turbines is possible. The next step is extract the utilization distribution (UD) from this track, effectively transforming the point data into a volume. As the length of the track is short (<48 h), the UD will be one of occurrence as opposed to a home range [62,64]. As the track is autocorrelated in time and space [65,66], the biased random bridge method implemented in the adehabitatHR package was applied [52]. Where data are taken to be representative of the individual's use of the space over a longer timespan (e.g., an entire breeding season), a home range estimator such as auto-correlated kernel density estimation may be applied instead, implemented in the 'ctmm' package [103,125].
Appendix A.5.1. Utilization Distribution The "BRB()" function in adehabitatHR accounts for autocorrelation by including a non-random "drift" component to the estimation of the probability density function [67]. The diffusion coefficient D (in units of the coordinates, in this case m 2 /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()", described in more detail in [69]. The minimum smoothing parameter hmin (m, the minimum uncertainty surrounding an animal's position) was set to 200 m following the distribution of distances between consecutive fixes [70]. Grid size will depend on study species' movement ecology (range of distances between fixes as well as maximum range), GPS interval and error, scale of environmental variables being compared with, and computing power. The smaller the grid size, the higher the resolution, the larger the file. In this case, grid size (100 m) was selected based on range of distances between fixes, scale of tidal lease site (<5 km), and reasonable file size. more detail in [69]. The minimum smoothing parameter hmin (m, the minimum uncertainty surrounding an animal's position) was set to 200 m following the distribution of distances between consecutive fixes [70]. Grid size will depend on study species' movement ecology (range of distances between fixes as well as maximum range), GPS interval and error, scale of environmental variables being compared with, and computing power. The smaller the grid size, the higher the resolution, the larger the file. In this case, grid size (100 m) was selected based on range of distances between fixes, scale of tidal lease site (<5 km), and reasonable file size. The 50% and 95% isopleths (i.e., contours) typically extracted from utilization distributions represent the 'core' and 'active' areas of use, respectively [64]. As we want to assess potential for overlap with turbine structures at the sea surface (e.g., floating mooring) and especially underwater (e.g., blades), known dive locations as determined by the TDR should be retained as far as possible, as these are the highest risk areas in terms of collision. Therefore, check whether 50 or 95% isopleths should be used in further analyses. This has ramifications for level selection in datasets where dive locations are unknown (e.g., GPS-only data). #get 50 & 95% isopleths vertex_50<-getverticeshr(ud, 50, unin =c("m"), unout=c("m2")) vertex_95<-getverticeshr(ud, 95, unin =c("m"), unout=c("m2")) Appendix A. 5
In the 50% isopleth, 0% of dive locations are retained, while virtually all are in the 95% isopleth. Therefore, the 95% isopleth should be used in this case. The utilization distribution is now ready to be used in further analyses (e.g., assessing overlap with marine renewable areas).

9, x FOR PEER REVIEW 25 of 32
In the 50% isopleth, 0% of dive locations are retained, while virtually all are in the 95% isopleth. Therefore, the 95% isopleth should be used in this case. The utilization distribution is now ready to be used in further analyses (e.g., assessing overlap with marine renewable areas).   In the 50% isopleth, 0% of dive locations are retained, while virtually all are in the 95% isopleth. Therefore, the 95% isopleth should be used in this case. The utilization distribution is now ready to be used in further analyses (e.g., assessing overlap with marine renewable areas).