A utilization distribution for the global population of Cape Vultures (Gyps coprotheres) to guide wind energy development

The rapid development of wind energy in southern Africa represents an additional threat to the already fragile populations of African vultures. The distribution of the Vulnerable Cape Vulture Gyps coprotheres overlaps considerably with wind energy development areas in South Africa, creating conflicts that can hinder both vulture conservation and sustainable energy development. To help address this conflict and aid in the safe placement of wind energy facilities, we map the utilization distribution (UD) of this species across its distributional range. Using tracking data from 68 Cape Vultures, collected over the last 20 years, we develop a spatially-explicit habitat use model to estimate the expected UDs around known colonies. Scaling the UDs by the number of vultures expected to use each of the colonies, we estimate the Cape Vulture population utilization distribution (PUD), and determine its exposure to wind farm impacts. To complement our results, we model the probability of a vulture flying within the rotor sweep area of a wind turbine throughout the species range, and use it to identify areas that are particularly prone to collisions. Overall, our estimated PUD correlates well with reporting rates of the species from the Southern African Bird Atlas Project, currently used to assess potential overlap between Cape Vultures and wind energy developments, but it adds important benefits, such as providing a spatial gradient of activity estimates over the entire species range. We illustrate the application of our maps by analysing the exposure of Cape Vultures in the Renewable Energy Development Zones (REDZs) in South Africa. This application is a scalable procedure that can be applied at different planning phases, from strategic, nation-wide planning to project-level assessments.


Introduction
With the global climate emergency, society needs to transition away from fossil fuels towards renewable energy sources. At the same time, there is a need to design strategies to minimise the environmental impacts associated with these alternative forms of energy production. Wind energy is one of the fastest growing forms of renewable energy (REN21 2021), but is known to result in the direct mortality of birds by collision with turbine blades, particularly of large soaring species (Thaxter et al. 2017;Drewitt and Langston 2006;Schuster, Bulling, and Köppel 2015).
Across Africa, wind energy offers opportunities to reduce national dependency on fossil fuels, increase energy security, and reduce greenhouse gas emissions. In South Africa, the aim is to establish 17.7 GW of installed capacity over the next decade, making wind energy one of the main alternatives to coal for the country (IRP 2019). South Africa also represents an important stronghold for the Cape Vulture Gyps coprotheres (Piper 2005); endemic to southern Africa, and with most of its population restricted to South Africa, Lesotho, and eastern Botswana (Piper 2005;Oikantswemang et al. 2021). Currently, the population of this Vulnerable species consists of 9600-12800 mature individuals, and has undergone severe declines since the 1990s (BirdLife International 2022). Although mitigating the effects of wind energy on Cape Vultures is considered of highest priority in South Africa, there have already been at least twenty fatalities to wind turbines in this country (Perold et al., 2020;Ralston Paton et al., 2017, unpublished data).
Understanding the spatial distribution of sensitive species is vital to ensure wind energy development is restricted within their core use areas. Tracking data have been used to develop predictive models for early identification of core areas for some raptors, and identify high collision risk locations (Murgatroyd, Bouten, and Amar 2021;Reid et al. 2015;Katzner et al. 2012;Poessel et al. 2018;May et al. 2013). These models can be particularly useful during early development planning to facilitate cost-effective risk assessments.
Over the last two decades, multiple organisations have attached GPS tracking devices to Cape Vultures with different objectives, such as: investigating habitat preferences, post-release survival of rehabilitated and captive bred individuals, or identifying causes of mortality (e.g. Bamford et al., 2007;Jobson et al., 2020;Kane et al., 2016;Phipps et al., 2017). In this study, we bring together most of these tracking data to estimate the utilization distribution (UD, Van Winkle 1975;J. Signer, Fieberg, and Avgar 2017) of the species across its entire range, to help wind energy stakeholders identify high-risk areas, and plan developments accordingly.
A UD measures utilization intensity and it can be used to estimate the probability of a bird using any region in space. Therefore, it can also be used to estimate the probability of encountering a wind turbine placed on any given region. We recognise that a UD alone does not inform collision risk, because it does not account for the probability of collision given an encounter with a wind turbine (Masden and Cook 2016;Murgatroyd, Bouten, and Amar 2021). Another limitation of a UD is that it provides information about the activity distribution of any one vulture, but it does not account for the number of vultures present in an area. To estimate the expected utilization distribution of the Cape Vulture population, we combine UDs around known colonies an roosts with information about the number of individuals using them.
More precisely, we use tracking data to build a behavioural model based on habitat preferences and movement patterns of Cape Vultures. Using this model, we simulate movement around breeding colonies and roosts, to map an estimated steady-state UD. We then scale this UD by the size of the local populations, as determined by the counts at breeding colonies and roosts.
Finally, we add all colony-wise scaled UDs to map the expected spatial distribution of the global Cape Vulture population at a given time.
Because of the recognised risk of wind energy collisions for this species, we are also interested in identifying areas where the species is likely to be flying specifically within the sweep of a wind turbine. However, because predicting flight height is challenging and adds an extra layer of uncertainty to utilization intensity estimates, we investigate UD both in two and three dimensions (2-D: UD at any height, and 3-D: UD at risk height). We demonstrate the applicability of our maps, in the context of strategic environmental impact assessment, by analysing the exposure of the Cape Vulture population within the Renewable Energy Development Zones (REDZ) in South Africa (DEFF 2019). This procedure can be adapted to the evaluation of wind energy developments at different scales (e.g., individual wind energy facilities), and to guide other forms of infrastructure (e.g., mining, airports, forestry, etc).

Modelling framework
We collated Cape Vulture tracking data from multiple sources to investigate the species' habitat preferences and ranging behaviour, and eventually estimate the distribution of activity of the population. There are several reasons why our raw tracking data are inadequate for directly calculating a utilization distribution, for example: i) birds were not tagged according to their numbers in the different colonies, and therefore might not be representative of the distribution of the population, and ii) birds with high resolution will have more weight on where tracking locations concentrate. Thus, instead of using raw data directly, we used it to model the general behaviour of the species, and then predicted its activity distribution throughout its range.
We used step-selection analysis to model space use (see Northrup et al., 2021), and logistic regression to model the probability of vultures flying at wind turbine height. In both models, we chose explanatory variables that were relatively in-variant over time to facilitate generalisable predictions. For example, we favoured topographical variables such as elevation, slope and terrain ruggedness, which are known to influence Cape Vultures movement, or soaring birds movement in general (Duriez et al. 2014;Hanssen, May, and Nygård 2020;Martens et al. 2020), over measures of orographic uplift (Shepard et al. 2013), which are more detailed, but harder to predict over time. In addition, the inclusion of climatic variables, such as temperature or wind speed to explain uplift, could require fitting complex variable interactions, and we were mindful of extra complexity in already complex models. We focused on prediction, rather than causal inference, and therefore model selection was based on predictive performance. Consequently, marginal effects of individual covariates should be interpreted cautiously.
We used the fitted models to simulate vulture movement around core colonies (see section 2.2 for the definition of core colonies), until a steady-state UD was reached. We estimated one UD for each core colony, and from these we calculated two distributions: i) the mean UD across individual colonies, which provides a measure of utilization intensity (or probability, if integrated over an area) across the distribution range of the species, disregarding colony size, and ii) the population utilization distribution (PUD), calculated by summing individual-colony UDs, after they have been scaled by the number of vultures estimated to use each colony (Figure 1). Note that we can switch between an unnormalized version of the PUD to obtain the number of individuals expected to use any location at any given time (population spatial distribution), and a normalized PUD that integrates to one over the whole Cape Vulture range to obtain a populationlevel UD (population utilization distribution), by dividing by the total number of vultures in the population. Both UD and PUD represent long-term averages over time.

Tracking data
Our tracking data came from 79 Cape Vultures tagged for different projects, and thus included multiple tag types with different fix intervals, ranging from a median relocation time of around one minute to nine hours. From these birds, we excluded those with less than one month of tracking duration, and those with the majority of fixes being stationary, based on the tag's speed readings (see identification of movement states in section 2.4). Irregular tracking sections in terms of resolution and/or geographical location (e.g., impossible, sudden changes in location) on the remaining birds, were visually detected and removed from the analysis. We removed fixes requiring the bird to travel >200 km/h on average. Some of the birds had been rehabilitated from injuries (see Appendix S1: Table S1). Although rehabilitated individuals could behave differently to wild birds (Jobson et al. 2020), we believe removing birds with short or irregular tracking periods, and those that spent too much time at their release site (which was the main concern of Jobson et al. 2020), excluded the most problematic individuals.
We calculated flight height by subtracting the altitude above sea level provided by the tags from the altitude at ground level provided by a 30 m resolution digital elevation model (NASA JPL 2013). After subtracting ground altitude, we obtained some negative heights, for two main reasons: i) for 13 of the birds, height readings systematically wrapped down to zero at a height of 2,042 m (birds flying at 2,042 m a.s.l. would appear to be at 0 m, those flying at 2,052 would appear to be at 10 m, and so on), and ii) by stochastic errors in tag readings.
To separate stochastic errors (typically small) from systematic errors produced by wrapping heights, we defined thresholds below which heights were probably wrapped. We did this by visually inspecting the height distribution for each bird (see Appendix S1: Figure S1). We then added 2,042 m to any height below this threshold. Other data processing specific to a particular modelling phase is described in the relevant sections below.
We classified tracked vultures as either adult, if it was older than four years, or non-adult. Nonadult individuals tracked into their fifth year were thus re-classified as adults.

Colony, roost and supplementary feeding site data
We collated information on the location and number of birds associated with colonies and roosts throughout the Cape Vulture distributional range. These data were collected by multiple individuals and organisations over an extended period (1970 -2020), with varying protocols, ranging from direct observation to helicopter surveys. Some of these surveys were published in a plethora of sources (see Appendix S1: Section S7), but many others appear only on official or private databases. Collating all these records was a massive undertaking by different authors of this manuscript over the years. To unify their records and form a common dataset was challenging, and we had to make subjective decisions, based on knowledge about the Cape Vulture population and colony monitoring experience. For example, we had to identify duplicated colonies with different names depending on the source, or colonies that were the same but located in slightly different locations. Those colonies that we were confident were active and correctly located, as well as roosts with average counts over time of at least 50 individuals, were used to as core colonies to define the main areas of vulture activity. Those with less certain data were used to locate areas that might be attractive, but not necessarily central, for the activity of the vultures (see section 2.3 to see how these colonies were used in the models).
To determine the expected size of a colony we took the average number of individuals counted in that colony over the years. Colonies typically had counts for the number of adults, breeding pairs or nests. To estimate the number of adults using the colonies, we doubled the number of breeding pairs or nests, when adult counts were missing. We made the distinction between adult and nonadult individuals using each colony, because they could have different behaviour (Martens et al. 2018). Thus, when the number of non-adults at a colony was missing, we estimated it following Piper (1993), who suggested that the proportion of breeding birds in the population is around 0.75 (we refer to breeders as adults and to non-breeders as non-adults). Additionally, we tested several population matrices based on survival and fecundity values presented in Borello & Borello, (2002), Piper, (1993, and Schabo et al., (2017), and the expected steady-state population, given by the main eigenvectors of the matrices (Caswell 2001). All suggested a roughly similar proportion of adults.
Tracked vultures were allocated to a specific point that represented its central colony, which was identified using the maximum kernel density estimate obtained from those fixes captured two hours before sunset and two hours after sunrise, when the vultures are most likely at the colonies (see Borello & Borello, 2002). Acknowledging that vultures (notably non-adults) could change their central colony, we repeated this process for each year of tracking.
The central Eastern Cape region presents high vulture activity levels based on field observations (Boshoff et al., 2009 and G.T. pers. comms), yet most of our records of colonies and roosts in the vicinity had low or no counts. To avoid completely missing activity in the area, for colonies with no counts, we provisionally allocated a colony size equal to the overall median colony size. To highlight the uncertainty associated with this area, we created a 50 km buffer (Venter, Martens, and Wolter 2019) around the convex hull enclosing the modified colonies that is presented in the final maps.
Supplementary feeding sites are locations where animal carcasses are placed for vultures to feed on. There are a number of these sites in South Africa that could influence Cape Vulture behaviour, and we incorporated these into our models. We used an updated data set of supplementary feeding sites based on the data by Brink et al. (2020).

Environmental data sources
We used a digital elevation model with a 30 m resolution (NASA JPL 2013) to obtain terrain elevation, slope, and a vector ruggedness measure to define ruggedness (Sappington, Longshore, and Thompson 2007). The distributional range of the Cape Vulture is topographically quite variable (e.g., high altitudes in Lesotho, low but rugged coastal regions in the Western Cape Province of South Africa or the rather flat terrain in Namibia). To account for different availability in topographic features, we divided the study area into four quadrants with the centre at (latitude, longitude): -27.5, 24.0. The quadrant identifier was used as a covariate in some of the models.
We used the 100 m Copernicus Land Cover data sets (CLC), spanning 2015-2019, to define land use (Buchhorn et al. 2020). Fixes during this period were annotated with land use data in the corresponding year, earlier fixes with data from 2015, and more recent fixes with data from 2019. The numerous land cover categories in CLC were reclassified to simplify the analyses, and avoid overfitting. Producing more general categories should also reduce the impact of annotating fixes with asynchronous land use data. Thus, all forests were classified as "closed" habitats, shrub and grass lands as "open" habitats, and cultivated land as "crops". Other categories include inland "water" and "urban" areas. To determine whether locations were within protected areas, we used the World Database of Protected Areas (UNEP-WCMC 2019).
All environmental variables were projected onto a ca. 1 km 2 grid (i.e., 0.01 degrees x 0.01 degrees), giving each grid cell the average value of the variable in the original scale, if the variable was numerical and the most abundant value if it was categorical. This same grid was used both for model fitting and predictions.

Step-selection model
We used a hierarchical step-selection model (Muff, Signer, and Fieberg 2020) to explore habitatselection and movement patterns. In these models, tracking fixes are compared with available locations (pseudo-absences). These pseudo-absences are sampled from the region considered available to the bird, which is conditional on observed locations, step-lengths and turning-angles (see Fieberg, Matthiopoulos, Hebblewhite, Boyce, & Frair, 2010).
We regularised (equally spaced fixes in time) each bird track to the highest possible resolution, and then modelled adjustments to the movement kernel following an integrated step-selection approach (Avgar et al. 2016;Fieberg et al. 2010;Forester, Im, and Rathouz 2009). Thus, the entire track of each bird was regularised to the same resolution, but different birds could be regularised to different resolutions. Consequently, pseudo-absence points for each bird were conditioned by the different step-length and turning-angle distributions. To explore systematic effects of resolution on selection coefficients, we also fitted a model that included interactions between habitat variables and tracking resolution.
Step-length was given a gamma distribution, and therefore, we tested step-length and the natural logarithm of step-length as covariates in the step-selection model to correct preliminary shape and rate parameter estimates (Avgar et al. 2016). Turning-angles were given a Von Mises distribution, so a cosine of turning-angle was tested as a covariate to correct preliminary concentration parameter estimates (Avgar et al. 2016).
Step-length and turning angle effects on step-selection could interact with tracking resolution to account for the fact that different time steps would have produced different preliminary parameter estimates for the gamma and Von Mises distributions. The likelihood of the step-selection process was approximated by a Poisson distribution (Muff, Signer, and Fieberg 2020).
We considered that distance from the central colony should be highly influential for Cape Vultures' ranging, because they are gregarious and central-place foragers. Thus, we tested as covariates both distance to central colony, to represent a linear effect of distance, and the natural logarithm of distance to the central colony, as an additional term that accounts for a short-range effect that changes rapidly with distance. Similarly, we predicted some potential attraction towards other colonies, and towards supplementary feeding sites (Schabo et al. 2017;Brink et al. 2020;Martens et al. 2020). We included an interaction effect between vulture age class and these distance covariates (to central colony, other colonies, and supplementary feeding sites), because non-adult birds tend to range further from colonies compared to adult birds, and behave differently around supplementary feeding sites (Monsarrat et al. 2013;Schabo et al. 2017). Sex was not included in the model because most tracked vultures were not sexed reliably. We also included terrain elevation, slope, and ruggedness, and their interaction with geographical quadrant, as well as land cover and protected areas. We incorporated random effects (i.e., random slopes) for individual vultures for the topographic covariates, as well as for those related to distance to colonies and supplementary feeding sites. We also included random intercepts for each step with a fixed large variance (Muff, Signer, and Fieberg 2020). We excluded random effects for other variables to reduce model complexity and avoid convergence problems.
We compared the predictive performance of different explanatory variables using a crossvalidation procedure. To assess each model, we divided the data into three parts, randomly allocating the same number of birds to each part, and using two parts for fitting the model (training set) and one part for validation (test set). We repeated this procedure ten times. For each bird in the test set, we took a sample of 500 steps with their associated pseudo-absences, to have all birds equally represented. After fitting the model to the training set, we predicted a selection score for the fixes and pseudo absences of the test set. Then, for each test-set bird, we ordered the data by increasing selection scores and divided them into ten equally-spaced bins. Our choice of using ten bins was arbitrary, but based on the fact that we would define ten activity levels in our final maps (see section 2.5). We then calculated the proportion of actual fixes to pseudo absences in each bin. If a model performed well, we would expect a higher proportion of fixes to pseudo-absences in those bins with high selection scores than in those with low selection scores.
We assessed model performance by looking at the Spearman rank-correlation between the bin rank based on selection scores, and the overall proportion of fixes to pseudo-absences associated with these bins ranks across all birds (Boyce et al. 2002).

Identifying movement states
For the flight height models, we used only fixes of moving vultures. To separate resting from moving states, we analysed the instantaneous speed readings provided by the GPS tags. We used Hidden Markov Models (HMM) with two states for birds that had resolutions of < 3-hour fix frequency. HMMs fitted to birds with coarser resolutions often resulted in convergence errors, which suggested there should be little temporal structure in these data. Therefore, we fitted uncorrelated mixture models to these birds' data. For both HMM, and mixtures, we used gamma distributions. We used time of day, and distance to central colony, to other colonies and to supplementary feeding sites, as covariates affecting the moving state. Time of day was fitted as a quadratic term, since vultures were expected to be more active in the middle of the day when uplift is usually greater.

Modelling flights at risk height
We modelled the probability of Cape Vultures flying at collision risk height, using altitude readings (telemetry tags' readings minus ground altitude, see section 2.2) from those fixes estimated to be in a moving state. We used a Bernoulli generalized linear model, with a logit link, and two possible outcomes: i) at-risk, when flight height <= 300 m, or ii) not at-risk, when flight height > 300 m. We used 300 m as the threshold because the highest turbine we know of is currently c. 250 m high (see https://www.ge.com/renewableenergy/wind-energy/offshorewind/haliade-x-offshore-turbine), and turbine sizes are increasing (see Therkildsen et al., 2021).
This higher threshold is conservative, and provides additional insurance against flight height errors from the GPS tags.
Height risk models included an autoregressive term to account for heights being more similar when occurring closer in time. Because of the different temporal resolution of the tags, we could not use a common autoregressive effect for all birds. Instead, we defined an effect that decays continuously with time as if the bird was flying at risk height at t-1, or otherwise.
Therefore, regularisation was not necessary, and we used the raw resolution of the tracks, unless fixes were < 15 minutes apart, in which case, we subsampled to this resolution. The choice of this threshold was arbitrary and based on our opinion that at higher resolutions, the autoregressive terms would undesirably dominate over other covariates, and would contribute little extra information.
We tested the predictive performance on risk height for the following variables: terrain elevation, slope, ruggedness, distance to colony, distance to supplementary feeding site, and the time of day (quadratic). We incorporated random effects for individual vultures for all covariates. Predictive performance was assessed using cross-validation, using a similar approach as previously applied to the step selection model. However, instead of using rank correlation for comparison, we predicted the probability of flying at risk height for each observation (conditional on the previous one) and calculated prediction performance using the area under the curve (AUC) of a receiver operating characteristic (Murtaugh 1996;Boyce et al. 2002). We could use AUC here because a fix is either at risk height or not, as opposed to in the habitat selection where pseudo-absences are not locations certainly unused by the vultures.

Simulating and mapping utilization distribution
Using our best step-selection model, we simulated vulture movement trajectories to create a UD around each core colony (see Figure 1). We defined a ca. 1 km 2 grid (i.e., 0.01 degrees x 0.01 degrees cells) across the species' range and simulated 2.4e 5 steps around each core colony. Using half this number (1.2e 5 ) produced near identical results, suggesting our simulations were sufficient to produce a fair approximation of the UD. Simulations started at each colony at sunrise, and produced hourly steps, with eight steps corresponding to one day. We chose this resolution to capture variation in activity due to time of day. Selection of consecutive steps was based on two probability distributions: i) a 2-D movement kernel defined by the step-length and turning angle distributions fitted to the birds that had one-hour tracking resolution (we used the average of the parameters across these individuals), and ii) a normalized habitat kernel produced by estimating habitat preference for each grid cell, based on the step-selection model and associated covariates. The product of movement and habitat kernels gave the resulting sampling probability for each cell. To account for variability in habitat preferences for different individuals, we sampled 40 coefficients from each random effect distribution (i.e., normal distribution centred at the estimated coefficient and variance equal to the estimated random effect variance), to represent 40 different birds. An equal proportion of the total number of steps allocated to each set of random coefficients.
Tracking data suggested that birds came back to their colonies regularly, particularly at night. To ensure we captured this pattern in our simulations, we set a probability of birds coming back to their colony that was dependent on how long they were away on a 'trip'. We defined a trip to have started when the bird was >5 km from the central colony and finished when it came back.
This distance was chosen arbitrarily, thinking that colonies spread over some area, and also to acknowledge that vultures make short distance movements while at the colonies. To estimate the distribution of trip durations, we fitted a negative binomial distribution to the number of days spent away from the colony (ndays) on each trip observed in the data. We used the probability of observing the number of simulated days away or less, at the end of each simulated day, as the probability of returning. We then sampled from a Bernoulli distribution with two outcomes: success, the bird came back; or failure, the bird stayed on the trip. Therefore, Where n and µ are the dispersion parameter and mean of the negative binomial, respectively.
Once we finished the simulations, we counted the total number of visits to each grid cell.
Additionally, we estimated the probability of each simulated step being at risk height using our flight-height model, and summed the probabilities of all simulated steps for each cell.
Considering each step being at risk height a Bernoulli trial with probability prh, the sum of the probabilities associated with each cell prh, provides the expected number of locations at-risk height.
The distribution of counts per grid cell obtained by simulation will still be noisy unless, an extremely large number of simulations are conducted, which is not practical due to the computational burden of our models. Instead, we smoothed the distributions fitting a Generalised Additive Model (GAM) to the number of visits to each cell (for 2-D) and to the expected number of locations at risk (3-D), using a Poisson likelihood. In this "smoothing" step, we included a similar set of covariates to those used in the step-selection model, but excluding those that had to do with time of day or movement (step-length and turning angles), because we were modelling the steady-state UD. Non-linear effects of minimum distance to any colony and supplementary feeding sites, tended to produce artifacts at the largest distances, so we used linear effects for these variables instead. We included non-linear effects for the covariates with only four knots to allow some flexibility, while retaining relatively smooth effects. We included a more flexible, penalised tensor product smooth between a cubic spline for distance to central colony, and a cyclic cubic spline for bearing from the central colony. This results in a spatial effect in polar coordinates taking the central colony as the origin. Finally, to produce a UD around each colony, we normalized the distribution of visits per grid cell to sum up to one over the entire range. The normalization ensures that the UD integrates to one, and it is a probability distribution. The previous steps produce a UD around each colony; to estimate the overall UD throughout the entire distribution range of the species, we calculated the mean UD across colonies (Figure 2). We used the mean to make the overall UD still be a probability distribution that integrates to one over the entire range.
We repeated the process above for adult and non-adult vultures, effectively, calculating two different UDs around each colony, and over the range of the species. To obtain a unified UD, we averaged adult and non-adult UDs.
In addition to estimating a UD covering the range of the species, we needed to account for different colonies being used by different numbers of vultures. To estimate the expected Cape Vulture population utilization distribution at any given time (PUD), we scaled each colony UD by the number of vultures at that colony. First, we multiplied UDs estimated for adults and nonadults, by the respective number of individuals, and then summed both scaled UDs to arrive at the distribution of the total population for each colony. We then summed all scaled UDs across colonies to estimate the global PUD. Additionally, we normalized the PUD so that it integrates to one, and represents population density (or proportion, if integrated over an area) rather than the actual number of vultures ( Figure 2). The result is effectively an average UD across colonies weighted by the number of vultures using each colony.
To delineate the final maps, we only used 99% of the accumulated PUD, because the remaining 1% represents vast areas where the presence of the vultures is only marginal.

Validating and applying the utilization distribution maps
As a crude validation, we compared our estimated PUD with Cape Vulture detection/nondetection data, at the pentad resolution (5' x 5' or ca. 9 km x 9 km) from the Southern African Bird Atlas Project 2 (SABAP2, M. Brooks & Ryan, 2020). We fitted a generalised linear model with a binomial error structure and a logit link function, to the SABAP2 reporting rates (number of detections in a pentad/number of visits, using number of visits as weights), and we used the mean cumulative PUD in each pentad as an explanatory variable. We used SABAP2 data collected between 2008 and 2020, and only those pentads with at least five visits. We used the cumulative PUD, because in early tests it correlated better than PUD or log(PUD) with SABAP2 reporting rates.
To illustrate the application of the PUD to assess encounter exposure, we analysed the encounter potential in the Renewable Energy Development Zones (REDZ), which are the South African government's preferential areas for renewable energy development (DEFF 2019). The procedure was as follows: 1. Define the broad reference area. This can be the entire distribution of the species, a country, or in general any region of interest. This reference will be used to give the assessment a context, in terms of the vulture population potentially affected. In our case, we defined both South Africa, for a national assessment, and the province the majority of each REDZs falls into, to provide a provincial assessment.
2. Define a polygon enclosing the locality of interest. In our study, each REDZ is a locality of interest. This could also be the footprint of a particular wind farm, for example, which could be useful for choosing the less risky among alternative sites.
3. Sum the PUD values of the pixels enclosed by the locality of interest, and by the reference region (if the assessment is to make that comparison).

Conduct relevant analyses, such as comparing Cape Vulture population exposure among
alternative sites or with the reference region.

Software implementation
We used the software R (R Core Team 2021) with the added functionality of the tidyverse The code used for processing and analysing the data (Cervantes 2023a) is available on ZivaHub, as well as a package (Cervantes 2023b) developed for simulating UD and PUD maps, and to update these maps as more information (e.g., updated colony counts) becomes available.

Tracking and colony data
From the 79 Cape Vultures tracked, 68 passed the first data-cleaning filters (29 adults and 39 non-adults, Appendix S1: Table S1), 13 of which were rehabilitated. Fixes span predominately across six countries: South Africa, Lesotho, Swaziland, Mozambique, Botswana, and Namibia ( Figure 3), between 2004 and 2020. We used 242,300 fixes for the step-selection analysis, and 287,030 fixes for the flight height analysis. The resolution of some birds had to be considerably reduced to achieve regular fix intervals for fitting step-selection models, because they had irregular sampling intervals (see Appendix S1: Table S1). For example, schemes with 2, 4 and 18 hour time intervals depending on time of day or night, resulted in regular time intervals of 24 hours (with 2 hour tolerance).
We identified a total of 729 breeding colonies and 220 roosting sites (Figure 3), of which 165 colonies and 5 roosts were considered core colonies. Colony size and roost size varied considerably, from colonies as small as 9 individuals to as large as 2,136 individuals (mean = 88.7, SD = 246). In total, there was a yearly average of 17,038 birds (of all ages) counted at all colonies combined. The mean proportion of adults to the total in our data was 0.76 (n = 24, median = 0.8, SD = 0.18).

Step selection model
We evaluated nine step-selection models (see Figure 4 and Appendix S1: Table S2). The crossvalidation Spearman correlation score of the best step-selection model was 0.979, suggesting very good performance at predicting habitat use (Boyce et al. 2002). The largest effects in this model ( Figure 5 and Appendix S1: Table S4) corresponded to distance from colonies, roost sites and supplementary feeding sites, suggesting vultures systematically moving towards these locations. However, the effect of these variables differed widely between individuals, as indicated by the large random effects' variance. In fact, we found possible avoidance of the central colony for some non-adult individuals. We found a positive selection of high elevations and ruggedness, and avoidance of urban areas and water bodies. The effect of elevation also differed markedly between individuals. Time of day interacted with distance to the central colony, as well as with step-length, suggesting longer displacements away from the central colony towards midday. Log of step length and cosine of turning angles had little predictive power, even at the highest (1-hour) resolution and were excluded from the final models.

Flight height model
We identified two states with the HMMs and mixture models. The "moving" state was clearly identified by the models and captured a range of speeds averaging 25-50 km/h (Appendix S2: Figure S2), depending on the bird. However, we noticed that the resting state, with averages of <5km/h might be capturing two different states: a pure resting state and a very slow moving state (Appendix S2: Figure S2). We considered that these short displacements over extended periods of time could be considered "resting" for our purposes.
The cross-validation score (AUC score) of the best flight height model, out of the six candidates, was 0.729, indicating an acceptable performance (see Hosmer, David W. et al., 2013, Figure 4 and Appendix S1: Table S3). In the best model, distance from the central colony, elevation, and slope had the largest effects (i.e., mean across individuals), with at-risk flights being more likely at higher elevations, on steeper slopes, and closer to the central colony. Effects of distance to colonies and supplementary feeding sites, as well as elevation, had the largest variation across individuals (largest random effects' variance).

Utilization and population distribution maps
Our predicted UD suggests that the largest utilization intensity is located in the Eastern Cape of South Africa, and south of Lesotho driven by the large concentration of colonies in these regions ( Figure 6). In addition, our PUD map, highlights the high population density around the large colonies in north-eastern South Africa and eastern Botswana. UD and PUD at rotor height were similar to the 2-D UD and PUD, but more concentrated around colonies (see Appendix S1: Figure S3).

Validation and application of the population distribution maps
SABAP2 reporting rates of Cape Vultures was highly positively correlated with our cumulative PUD estimates (Coef. = -7.59, SE = 0.061, p < 0.001). Note that PUD accumulates from areas of high utilization intensity to areas of low intensity, and thus the inverse relationship with reporting rates. Locations where Cape Vultures were more likely to be observed at SABAP2 pentads, were also predicted by our models to have higher vulture utilization (Figure 7 and Figure 8).
Examining the population utilization distribution accumulated by the eleven South African REDZ, we found that only three of them overlap substantially with the Cape Vulture population ( Figure 9 and Appendix S1: Table S6

Discussion
Using tracking data from 68 vultures, and the location of 167 core colonies, we estimated and mapped the expected UD of the Cape Vulture across its distribution range. Complemented with colony counts, the UD allowed us to estimate the exposure to wind turbines of the Cape Vulture population. Our predictive habitat use model, on which the UD was based, performed well based on cross-validation rank correlation with left-out data. Additionally, we incorporated predictions from a model for the probability of flights at risk height, with an acceptable predictive performance, although not quite satisfactory.
Considering the costs of wrong predictions, and because additional modelling incorporates more assumptions and potential errors, we prefer a more conservative approach and favour our 2-D distribution maps over the 3-D maps, for initial site screening and strategic impact assessments.
The opportunistic foraging strategy of the species may lead to unexpected flights at low heights (<300 m), even in areas predicted to be relatively safe in terms of flight height, if the right conditions occur (e.g., presence of carrion). Although these situations might be infrequent, the consequences can be fatal, particularly considering that multiple individuals could be involved in these events, given the gregarious nature of this species. However, used with caution, the 3-D maps may still be useful for cases in which an extra level of detail is desired. For clarity and unless otherwise stated, we focus on the 2-D version of the UD/PUD.
Our modelling approach favoured predictive performance rather than causal inference; and therefore, specific effects of individual variables on UD must be interpreted cautiously. Keeping this in mind, and to help with the interpretation of our distribution maps, we note that according to our step-selection model, the main driver of vulture activity is the location of colonies, roosting sites and supplementary feeding sites. Similar patterns of association of Cape Vulture movements with colonies and supplementary feeding sites were described by Kane et al., (2016), with the caveat that some of their data were also used in the present study. They described some variability in preference for areas close to colonies and supplementary feeding sites across individuals, which is also captured by (the random effects of) our model. The simulated UD incorporates variability across individuals using a sample of random coefficients; however, other patterns discussed by Kane et al., (2016) and Pfeiffer et al., (2015), such as seasonal variation in home range size, are not incorporated in our maps because we focussed on the average distribution across all seasons to guide general wind energy planning, rather than aiming to create season-or condition-specific dynamic maps.
Although time of day is a dynamic variable that is averaged out in our model, it is important to note that differential step-lengths throughout the day also have an impact on the emerging longterm distribution of activity around colonies and roosts. According to our model, step-length, which is indicative of displacement, tends to be shorter during the early and late hours of the day, when vultures are often found at colonies and roosts (see Borello & Borello, 2002). Shorter displacement during these hours results in the accumulation of activity around colonies and roosts. High-risk flights (under 300 meters) also tend to occur more frequently around breeding colonies, which is to be expected, given that vultures tend to land and embark on social interactions at these locations.
Currently, the guidelines for impact assessments of wind farms in relation to Cape Vultures (Pfeiffer and Ralston-Paton 2018), recommend SABAP2 data as a screening tool, distinguishing between areas with reporting rates <14% (low-density), and areas with reporting rates >14% (high-density). Other areas close to colonies, roosts or supplementary feeding sites, together with topographical features that could increase collision risk, are also classified as high-density areas.
Proposed wind energy projects in high-density areas require more intensive monitoring prior to approval. Our tool offers an alternative screening tool, combining all of this information into a common modelling framework, and at a finer resolution (ca. 1 km vs ca. 9 km). In addition, our maps provide a gradient of activity intensity, and a spatially-explicit abundance approximation, rather than just occurrence, and cover the whole range of the species, as opposed to SABAP2, which has uneven coverage, as well as some gaps in poorly surveyed areas. Coverage could admittedly be dealt with using other modelling techniques such as occupancy modelling (Altwegg and Nichols 2019), yet this would still result in occupancy probabilities.
Although, in general, there is a positive relationship between our estimated PUD and SABAP2 reporting rates, we see that the simple logistic regression used for the comparison leaves a lot of unexplained variation in reporting rates, and some discrepancies between PUD and reporting rates ( Figure 8). These disagreements could be explained by the non-linear relationship between abundance and detection rates (Royle and Nichols 2003;Lee, Fleming, and Wright 2018), and possibly poor coverage of certain areas by SABAP2, which could also decorrelate reporting rates from the underlying vulture activity. On the other hand, our risk estimates are model-based, and therefore subject to our assumptions that, if unrealistic, could also produce local disagreement between the estimated PUD and the actual vulture activity. Our recommendation would be to consider both sources, SABAP2 and our PUD maps, as complementary rather than mutually exclusive. Areas with high reporting rates and areas with high estimated PUD should both be examined carefully before planning wind energy development.
Whilst we believe that we have created a useful screening tool in the context of wind energy development for Cape Vultures, we also recognize that our maps have limitations. The UD generated by our models is influenced mainly by the location of core colonies, and the PUD also by the size of these colonies. Therefore, our predictions are highly sensitive to the accuracy of our colony data. Whilst we attempted to collate all known data on colonies and roosts, we used data from multiple sources, which led us to make subjective choices when integrating them into a unified database. We also focused on the average number of individuals counted at the colonies over the years, which may not be representative of the current situation. In addition, we identified certain areas within the species range where data are less comprehensive. For example, we highlight that our predictions from the central Eastern Cape are less reliable, based on the mismatch between our colony data, and field observations, which recorded higher activity than expected if these data were accurate (Boshoff et al., 2009 and G.T. pers. comms). We find two possible explanations for this: i) the area is heavily used by non-breeding individuals (Boshoff et al. 2009), which results in small breeding colonies, but seasonal, and possibly extensive, roosting in the area, and ii) these sites are more difficult to survey because they are smaller and more widespread than in other areas.
Another limitation that stems from the central role of colony data in our model, is that all our records correspond to locations at cliffs. Although predominantly a cliff roosting species, Cape Vultures do occasionally roost in trees and on powerline pylons in parts of their range (Bamford et al. 2007;Phipps et al. 2013). Ignoring this plasticity, may result in missing unusual roosts, potentially dispersed over large areas. Although these should represent a relatively small fraction of roosts, with the current model structure, UD/PUD maps could be improved by providing, if known, the centroid, and the total count of areas with important concentrations of Cape Vultures that are not associated with cliff roosting sites.
Finally, vulture counts at colonies focus primarily on counting adults and breeding pairs. The number of non-adults present at the colonies often had to be derived from adult counts, based on a simple ratio (3:1 adult to non-adult). This ratio is potentially inaccurate, although we observed a similar average ratio across colonies on our limited data. The dispersion observed in our data also suggests that the proportion of adult to non-adult birds is not a constant, and varies widely between colonies. Although it is possible that this proportion depends on the size and location of the colony (Piper 1993), this was not clear from our data, and further research could help clarify if this type of pattern could be established. Nevertheless, some basic tests (see Appendix S1: Section S8) show that the distribution patterns emerging in our maps are quite insensitive to the age ratio, and to some other rather arbitrary choices we made during the simulations (e.g., minimum distance from the central colony that determines the beginning of a trip).
To have a sense of how accurate our colony counts were, we summed the number of individuals in all colonies to estimate the total population of Cape Vultures. We obtained approximately 17,000 birds, 11,300 of which are mature individuals, which is a plausible number, within the range of the previous estimate of 9600-12800 mature individuals (BirdLife International 2022).
Regardless, we are hopeful for better data for some of these colonies in the future, and perhaps to some degree of standardization in colony counting. We are also aware of some ongoing initiatives to improve the coverage and accuracy of counts in the Eastern Cape of South Africa, where our uncertainty was greatest (GT pers. comm.). With that expectation for better data, we have developed the R package vultRmap (Cervantes 2023b) that can update our distribution maps, as new colony or roost information becomes available.
We demonstrate the application of our PUD maps analysing the current exposure of the Cape Vulture population within the designated REDZs in South Africa. Our analysis revealed that three REDZs (Stormberg Wind, Emalahleni and Overberg Wind) overlap substantially more than any others with Cape Vulture distribution, and together expose a potentially significant portion of both the provincial and global population to wind energy developments. Unfortunately, what represents a significant portion is not well defined currently; research on population viability of the species would greatly enhance the utility of these assessments. For now, we conclude that the Stormberg REDZ on its own could expose as much as 1.6 % of the global population, and as much as 6.7 % of the provincial population of the species. A relatively large percentage of the provincial populations are exposed in the other two REDZs (13.7 % and 12.3 % for Emalahleni and Overberg Wind, respectively). Interestingly, the PUD overlap with the Overberg Wind is concentrated in specific areas within the REDZ, as opposed to Emalahleni, where the exposure was more dispersed (Figure 6). This situation suggests different management strategies. The larger variance in the PUD at Overberg Wind suggests that avoiding specific areas with the highest exposure, could significantly reduce the potential risk of this REDZ for the Cape Vulture population. This is in contrast with Emalahleni, which presents a narrower range of PUD values; thus, making it more difficult to avoid potential risk, by careful siting of wind energy facilities.
We note that some of the REDZs intersect with several provinces, and that we have allocated them to the province covered the most. Although ecologically it might not make sense to use administrative boundaries, we still believe that this application illustrates the intended use of our PUD maps.
This illustrative procedure can also be applied to evaluate encounter potential at different scales or for different types of projects. For example, using an individual colony UD, we could assess the probability of encountering any one individual from that specific colony; if we multiplied this probability by the number of vultures using the area, we would arrive at the expected number of vultures encountered (or exposed). In contrast, using the range-wide mean UD, we assess the probability of encountering any one individual vulture across colonies, and using the UD scaled by the number of vultures at the colonies (i.e., PUD) we would assess the population exposed to encounters. The main difference between UD and PUD, is that UD does not incorporate information about the size of the colony, and therefore the utilization distribution is mostly influenced by the landscape, whereas PUD is also directly affected by the size of the colonies, with larger colonies concentrating higher utilization/exposure.
We believe that the set of maps we have produced can be a valuable tool for planning not only wind energy developments, but also other forms of infrastructure and land use, such as the delineation of protected areas or vulture safe zones. We also hope that our work complements other quantitative analyses, such as population viability assessments of the species, by providing a spatially-explicit estimate of the number of vultures potentially affected by human infrastructure. Lastly, we would like to underline the open-ended nature of this project and encourage further collaboration to improve these maps and build a more sustainable future for the Cape Vulture. We recommend that our approach could also be applied to other colonial vulture populations occupying areas where wind energy development is planned, for example for Rüppell's Vulture (Gyps rueppelli) in East Africa.     For the step-selection model estimates are separated for adults (ad) and non-adults (juv). The dots represent coefficient estimates, thick bars cover a 95% confidence interval and the thin vars cover a 95% of the distribution of the random effects (centered at the estimate). sl_ -step length, dist_col -distance to central colony, log_dist_col -natural log of distance to central colony, dist_col_any -distance to any other colony, dist_sfs -distance to supplementary feeding sites, elev -elevation, slope -terrain slope, rugg -ruggedness, prot_area -protected area, closedclosed habitats, crops -cultivated land, water -land covered by water, urban -urban areas, ttnoon -time to noon (10 am GMT), ttnoon_sq -time to noon squared, phi -autoregressive term. Figure 6. Two-dimensional utilization distribution (UD) and population (utilization) distribution (PUD) throughout the distribution range of the Cape Vulture. The colour gradient represents cumulative distribution levels, such that lighter colours enclose a larger proportion of the distribution than darker colours. The greyshaded polygon on the Eastern Cape province of South Africa encloses those colonies that were allocated a size equal to the median colony size for the country, due to lack of reliable colony size data for all colonies in this region. Within this region our UD and PUD estimates are less reliable. Note that it is the cumulative distribution of activity that is displayed in the maps (e.g., the contour around the 0.1 area encloses 10% of the activity). Figure 7. Comparison of the distribution of SABAP2 reporting rates per pentad (right) and mean cumulative population utilization distribution (cumPUD) per pentad (left). The grey areas represent those pentads with less than five SABAP visits, and were not used for the analysis.

Figure captions
Note that PUD accumulates from areas of high intensity to areas of low intensity; therefore, areas with cumPUD close to 1 have less utilization intensity than those with cumPUD close to 0. prediction interval. Note that PUD accumulates from areas of high intensity to areas of low intensity; therefore, areas with cumPUD close to 1 have less utilization intensity than those with cumPUD close to 0. Figure 9. Cape Vulture population utilization distribution (PUD) accumulated by (exposed at) South African REDZs. In the top panels, violin plots show the distribution of log(PUD) values within each REDZ (left) and average PUD value ± standard deviation (right), to aid interpretation. Note that the minimum of the bars has been truncated at zero to emphasize that PUD cannot be negative. In the lower panel, the location of REDZs overlaid with the estimated PUD. Note that it is the cumulative distribution that is displayed in the maps (e.g., the contour around the 0.1 area encloses 10% of the activity).