Bird-window collisions are a leading cause of direct anthropogenic avian mortality, yet our state of knowledge regarding this threat relies heavily on eastern North American studies. Seasonal patterns of collision mortality may differ along the Pacific coast, and western North American species remain understudied. We therefore surveyed a stratified random sample of 8 buildings for collisions at the University of British Columbia, Vancouver, Canada over 45-day periods during 2 winters, 1 spring, 1 summer and 1 fall season between January 22, 2015 and March 15, 2017. After accounting for the rate of scavenging and efficiency of observers in finding carcasses, we estimated that 360 collision fatalities (95% C.I.: 281 to 486) occurred over 225 days of monitoring. Collision mortality was highest in fall, but in contrast to most published research, collision mortality was intermediate in both winter and spring, and was lowest in summer. In winter 2017, we performed point count surveys to assess whether individual species are disproportionately vulnerable to collisions when accounting for population size, and found that the Varied Thrush (Ixoreus naevius) was 76.9 times more likely to collide with buildings, relative to average species vulnerability in winter. To our knowledge, this is the first study to report the Varied Thrush as a species that is disproportionately vulnerable to collisions. Further studies are needed to assess the vulnerability of Western North American species and subspecies and to determine whether similar patterns of seasonal collision mortality are found elsewhere.

Study Area

We conducted our research at the 420-ha Vancouver campus of the University of British Columbia (hereafter referred to as UBC), along the Pacific coast of southern British Columbia (BC), Canada (49.2606° N, 49.2606° W). The campus is situated within the Fraser Lowlands (elevation 0 - 152 m) bordered by the Strait of Georgia and the Coast Mountains (Figure 1). UBC is located on a peninsula, separated from the City of Vancouver by Pacific Sprit Park, a 874-ha park comprised mainly of 70 - 130 year-old second-growth forest (Metro Vancouver 2019). Tree canopy at UBC is predominantly coniferous; dominated by western red cedar (Thuja plicata), Douglas-fir (Pseudotsuga menziesii), western hemlock (Tsuga heterophylla), and grand fir (Abies grandis). ~58% of campus is in park-like open condition, with mowed grass, a number of coniferous and evergreen broad-leaf plantings, and scattered deciduous trees; primarily red oak (Quercus ruba) and maple species (Acer spp.) (Dyck 2016). The campus and surrounding forest lies within the Coastal Western Hemlock (CWH) Biogeoclimatic Zone (DataBC 2018) characterized by warm dry summers and mild wet winters, with a climax forest dominated by western hemlock (Tsuga heterophylla), as well as Douglas-fir (Pseudotsuga menziesii) and western red cedar (Thuja plicata) (Pojar et al. 1987, 1991). Broadly, the region is part of the temperate rainforest biome, also known as the Northern Pacific Rainforest Bird Conservation Region, composed of similar climate, habitat types and bird species, extending from the Gulf of Alaska, through to Northern California (Alaback 1996, North American Bird Conservation Initiative 2000, Rich et al. 2004).

Building Selection

We chose 8 buildings for bird-window collision monitoring at UBC, using a randomized design, stratified by building height and extent of surrounding vegetation (Figure 2; see Hager and Cosentino 2014, Hager et al. 2017). We used a map layer of 229 UBC building footprints and size-related attributes of all buildings >1 story, excluding parking garages (University of British Columbia 2015), to divide buildings into two categories: “short” buildings (2-4 stories) and “tall” buildings (5-12 stories). We used ArcMap 10.3.1 to generate a 50-m buffer around buildings and digitized all vegetation using the World Imagery base map (ESRI 2015) to calculate % vegetation within 50 m of each building. Within each building height class, we then randomly selected 2 buildings with >50% vegetation cover within 50 m of the building and 2 buildings with <50% vegetation cover within 50 m, as defined in Hager and Cosentino (2014). Excluding 1-story buildings eliminated most maintenance sheds and utility buildings with few windows. We further specified that study buildings must be a minimum of 100m apart, which resulted in one building being removed from the initial selection. We excluded a second building from our initial selection due to construction fencing which prevented access to the majority of the building perimeter.

Collision Surveys

We followed protocols outlined in Hager and Cosentino (2014) and Hager et al. (2017) to conduct a total of 155 bird-window collision surveys between January 22, 2015 and March 15, 2017 at our 8 study buildings. Surveys were conducted over 5 45-day periods consisting of 2 winter seasons (January 22, 2015 - March 7, 2015 and January 23, 2017 – March 15, 2017), 1 fall season (September 9, 2015 - October 23, 2015), 1 spring season (April 15, 2016 - May 29, 2016 and 1 summer season (June 16, 2016 - July 30, 2016). We aligned our spring and fall collision monitoring with the calendar dates when a local banding station records the highest numbers of migrants (WildResearch 2015). We performed a clean-up survey the day before each monitoring period to remove all evidence of bird-window collisions that accumulated prior to survey periods, including carcasses and feathers. In winter 2015 we surveyed every second day for 45 days. For the second winter season and all other seasons we conducted daily surveys for 21 days followed by surveys every second day for an additional 24 days; for a total of 45 days. We stopped surveys for 6 days during the second winter season due to a heavy snowfall that prevented us from reliably finding collision evidence. Surveys resumed after a clean-up survey to remove any evidence that accumulated while surveys were halted, and the survey period was extended accordingly.

We conducted collision surveys in mid-afternoon as recommended by Hager and Cosentino (2014), since window collisions often peak in the morning and continue at a lower rate later in the day (Klem 1989, Kahle et al. 2016). Two surveyors worked concurrently, walking around the perimeter of buildings in opposite directions, searching for carcasses within 2 m of the building façades. We defined a façade as the exterior face of a building, generally facing the same direction. In some cases, the complexity of a building footprint required us to consider several portions of a building footprint as one façade, e.g., several sides of an alcove, if we could not reliably discern which portion of the façade that a bird may have hit. One surveyor, ANP, remained constant throughout all 5 collision monitoring periods, while the second surveyor(s) differed seasonally. Evidence of a collision included whole carcasses and partial carcasses (e.g. wing, feet, bill and/or feather piles) (Hager and Cosentino 2014, Hager et al. 2017). We conservatively specified that feather piles must include a minimum of 10 tail, wing or body feathers confined to 50 cm diameter circular area (Ponce et al. 2010). Only two stunned/injured birds were found during surveys, and both died shortly after discovery. We removed all carcasses immediately following each survey to ensure double counting did not occur, and identified all intact carcasses to species. We also collected partial carcasses, including all feathers, and where possible these were later identified to species-, genus- or family-level, using online resources (U.S. Fish and Wildlife Service Forensics Lab 2018), hard copy guides (Scott and McFarland 2010), or by comparing to specimens from the UBC Beaty Biodiversity Museum.

Correcting For Biases Due To Scavenging And Missed Carcasses By Observers

Collision frequency estimates can be negatively biased due to carcass removal by humans and scavengers prior to their observation by surveyors, and imperfect detection of carcasses by surveyors (Klem et al. 2004, Bayne et al. 2012, Riding and Loss 2018). Therefore, it is important to design trials to estimate bias parameters when comparing building mortalities among studies, or among seasons and regions of interest (Loss et al. 2015). To account for carcass removal and imperfect detection, we conducted carcass persistence and searcher efficiency trials, then incorporated the results into statistical models using the Generalized Mortality Estimator in R package GenEst (Dalthorp et al. 2018a, b) to correct our collision mortality estimates (details below). We tested both intercept-only and season-specific models, but had insufficient data for both carcass persistence and searcher efficiency trials to allow us to apply bias corrections at the species-level, by building perimeter substrate or by other potential covariates.

Carcass persistence (CP). We conducted 103 carcass persistence (CP) trials to quantify bias due to scavenging of carcasses. Trial periods were chosen based on field worker availability and logistics of conducting frequent carcass checks; February 16 - March 7 in winter 2017 (n = 33), April 20 - April 28 in spring 2016 (n = 28), July 12 - July 30 in summer 2016 (n = 28), and September 10 - October 16 in fall 2015 (n = 14). We used carcasses from a range of previously window-killed birds, ranging in size from 4g (Anna’s Hummingbird, Calypte anna) to 128g (Steller’s Jay, Cyanocitta stelleri). Carcasses were stored in deep freeze and held under our Canadian Wildlife Service Salvage Permit, as outlined in the Ethics Statement. We placed thawed carcasses at a range of times throughout the day at a randomly selected study building and at a randomly-selected façade. The location along each façade was also chosen randomly, but was no closer than 2m from a corner or edge of an adjacent façade. Carcasses were placed on a variety of substrates, depending on the randomized location chosen, including concrete, asphalt, river rock, bark mulch, grass, or under bushes and ground cover. Collision monitors were alerted to the carcass persistence study and to the location of CP trial carcasses. However, as a precaution, data were checked carefully following all collision surveys to ensure that carcasses from CP trials were not inadvertently included in the collision mortality dataset. We checked carcasses every 5 – 6 hours after initial placement on the first day until early evening, in early morning and mid-day on the second day, and at mid-day on subsequent days until carcasses were no longer present (presumably removed by scavengers or campus maintenance staff), or until the end of the study period, whichever arose first. Carcass removal time was the interval between placement of a carcass at a façade until it was removed by a “scavenger” and no longer detectable by observers, as defined by Riding and Loss (2018).

We used bias-correction functions within the R package GenEst to fit a carcass persistence (or parametric survival) model to estimate the amount of time a carcass would persist, given the conditions under which it arrived (Dalthorp et al. 2018a, b). We tested whether a persistence model that explicitly varies depending on the season, or one that depends only on the intercept, best described the detection probability of a carcass over time. We fit one-parameter exponential and two-parameter (location and scale) Weibull models, and tested whether persistence parameters varied categorically by season. The location parameter in the exponential and Weibull models allowed covariates to shift the mean of the probability density function, whereas the additional scale parameter in the Weibull model allowed covariates to affect the degree of spread of the distribution. Since the sample size was small (n/K < 40) we chose the model with the lowest 2^nd-order AIC (AICc) goodness-of-fit statistic across models.

Searcher efficiency (SE). We estimated searcher efficiency by having a non-surveyor place a total of 29 thawed bird carcasses during the winter 2015 (n = 1), winter 2017 (n = 20), spring 2016 (n = 3), summer 2016 (n = 3), and fall 2015 (n = 2) study periods. Carcasses ranging in size from 4 g (Anna’s Hummingbird) to 160 g (Northern Flicker, Colaptes auratus) were placed at randomly-selected buildings, façades, and façade locations as described above for CP trials. Carcasses were removed from the trial if they were detected by at least one of the surveyors. Carcasses that were not detected on the first day were checked prior to each subsequent survey until they were located by observers, removed by a scavenger, or until the end of the collision monitoring period, whichever arose first. We used the R package GenEst to model searcher efficiency as a function of the probability that carcasses are located by observers on successive searches (Dalthorp et al. 2018a, b). The probability of finding a carcass on the first search after carcass arrival was denoted by the parameter, p, and k was the change in probability of finding a given carcass on the second and subsequent searches until the end of the collision monitoring period (Dalthorp et al. 2018a). We fit two models to estimate the observer bias parameters, including a model that depended only on the respective intercepts, and another model that allowed searcher efficiency to vary by season. We selected the best model by comparing the goodness of fit of models using AICc.

Bias-corrected mortality estimation: a two-stage parametric bootstrap approach. We used the function estM in the R package GenEst to produce a bias-corrected estimate of mortality (Dalthorp et al. 2018b). EstM estimated each carcass’ contribution to the total mortality in each search interval by accounting for the carcasses missed due to scavenging/imperfect carcass persistence (CP), and imperfect detection/searcher efficiency (SE) (Dalthorp et al. 2018a, b). We assumed that most birds that collided with buildings fell within 2 m of building façades (Hager et al. 2013, Hager and Consentino 2014, Riding et al. 2020), and therefore set the parameter density weighted proportion (dwp) or the proportion of carcasses expected to arrive in the area searched by observers (Simonis et al. 2018), to 1. The model uncertainty from the two bias correction models (carcass persistence and searcher efficiency models) was propagated into estimates of total mortality using parametric bootstrap methods (Efron and Tibshirani 1986). Specifically, function estM in the R package GenEst drew 1000 sets of estimated carcass counts from the asymptotic distributions of the maximum likelihood estimators of parameters from the two bias-correction models. The result was a matrix of carcass contributions where each row corresponded to a carcass, and each column corresponded to a simulated realisation of that carcass’ estimated contribution to total mortality. These simulated outputs resulted in bias-adjusted mortality estimates that included corrections for differences in search interval (time-increments between searches on the same building façade), corrections for seasonal difference in carcass persistence, and a correction for imperfect searcher efficiency.

Assessment of Relative Vulnerability to Collisions in Winter 2017

After 4 seasons of collision monitoring, we noted that some species were particularly prevalent in the collision mortality data. However, assessment of relative vulnerability to collisions requires an estimate of population size of species within the study area, as well as species-specific mortality counts. Therefore, in winter 2017, we conducted avian point count surveys during our collision monitoring period, and corrected these counts for imperfect detection. We used raw collision mortality counts to calculate relative vulnerability of species to window collision mortality on UBC campus as we did not have sufficient data to generate species-specific bias-corrected mortality counts as described above. To determine how population size influenced the likelihood of collision by species, we applied a two-step approach: 1) we estimated the change in total number of birds by species over time in response to population-level processes of mortality, immigration and emigration, and then; 2). we regressed mortality estimates on the corrected population counts to determine a species-specific vulnerability rating.

Winter 2017 point count surveys. We conducted 5-minute point counts at a randomly chosen façade of each building approximately weekly between January 25, 2017 and March 12, 2017, for a total of 72 point counts over 9 survey days throughout the winter 2017 collision monitoring period. We minimized factors that could contribute to detection error by: 1) recording all birds seen or heard by call notes within a small area to maximize the probability of detecting all individuals (i.e., a 50-m radius semi-circle, from the edge of building façades), and excluding birds that were only seen flying over; 2) conducting surveys during the morning hours when human foot traffic was low, and only in weather conditions favorable to maximum detection of all individuals (i.e., with no precipitation and winds below 13 km/hr), and; 3) using the same observer (KLD) to conduct all point count surveys to minimize variation in observer bias.

Population size estimates. We calculated the population size for each species using the Dail–Madsen model (Dail and Madsen 2011), which is a generalized form of the Royle (2004) binomial N-mixture model developed for open populations (i.e., the closure assumption is relaxed, allowing abundance to change over time). We chose this model for our winter point count data because it 1) assumes that abundance patterns are determined by an initial territory establishment process followed by gains and losses resulting from mortality, movements into the population (recruitment) and movements out of the population, and; 2) accounts for imperfect detection probability. For all species, we assumed that individuals could move across all sites (building locations) and all days within the winter 2017 collision monitoring and point count period. Thus, the number of gains across days was directly proportional to the number of survivors (individuals occupying winter territories), so we applied an auto-regressive model to recruitment, γ*N_s[i,t-1]; where γ = recruitment, N = the number of individuals per site, i, for each species, s, on survey day, t. As a result of the relatively short detection distance, we used a starting value of 1 for probability of detection. We applied a Poisson distribution to all species in which the total number of individuals detected across all 9 surveys, N_stot, was greater than 10. For species with N_stot < 10, we applied either a negative-binomial or zero-inflated Poisson distribution to a simple intercept-only model, and restricted immigration to zero. We then used the parameter estimates from this simple model as starting values for the more complex, open count model that allowed immigration across survey locations and days. For all models, we set the upper bound for discrete integration, K, at the greater value of N_stot + 1, or 10. We summed across all sites the resulting site-specific population size estimates corrected for detection probability and movement into the study area, divided by the number of survey-days (9), and multiplied by number of buildings surveyed to obtain population estimates for each species, N_corr. We calculated the root mean squared error on the residuals of each model to estimate residual standard error for each species’ estimate of N_corr. We estimated all parameters and obtained corrected population size estimates using the function, ‘pcountOpen’ and back-transformed the estimates for detection probability from each model using the function, ‘plogis’ in the R package, unmarked (Dail and Madsen 2011, Fiske and Chandler 2011, Hostetler and Chandler 2015, R Development Core Team 2018).

Calculation of relative vulnerability to collisions by species. Next we used our raw mortality counts per species for winter 2017 and species population estimates at our 8 study buildings (N_corr) to assess relative vulnerability of species to collisions during the winter 2017 season. We fit a simple linear regression to log₁₀ (X + 1) of the raw winter 2017 mortality counts on log10 of the population size estimates, N_corr (Arnold and Zink 2011, Loss et al. 2014). The use of the (X + 1) transformation on mortality counts allowed us to account for zero mortality for species that were observed in point counts but not in the mortality counts. We set regression coefficients for population size, N_corr to 1.0, calculated residuals, then raised 10 to the power of the absolute value of the residuals, for all species observed in our point count surveys. A slope of 1.0 assumes that mortality due to collisions for a given species is proportionate to their population size, as described by Arnold and Zink (2011). Vulnerability values represent the factor by which a given species is more (+) or less (–) vulnerable to collisions relative to the expected value for a species with average collision risk (residual = 0).