Risk of lethal vessel strikes to humpback and fin whales off the west coast of Vancouver Island, Canada

Vessel strikes are a source of mortality and injury for baleen whales, which can have population-level impacts. Spatial analysis of whale and marine traffic distributions provides a valuable approach for identifying zones with high collision risk. We conducted 34 systematic aerial surveys to estimate humpback Megaptera novaeangliae and fin whale Balaenoptera physalus densities off the west coast of Vancouver Island, Canada, including approaches to major shipping lanes in Juan de Fuca Strait, a gateway to the ports of southern British Columbia and Washington State. To predict whale densities, we fit negative binomial generalized additive models (GAMs) to sightings data, incorporating survey effort as an offset, and depth, slope, and latitude as environmental covariates. Humpbacks were primarily observed on the continental shelf, with highest predicted densities along the shelf edge (~200 m isobath), whereas fin whales were primarily distributed west of the shelf break (>450 m depth). We combined GAM-predicted whale densities with vessel traffic data to estimate the relative risk of ship strikes. Since vessel speed is an im portant determinant of lethality, we also calculated the relative risk of lethal injuries, given the probability that a collision occurs. Humpbacks were most likely to be struck along the shelf edge, the inshore approaches to Juan de Fuca Strait, and within the strait itself. Fin whales were most likely to be struck in the offshore approaches to Juan de Fuca and inside the western portion of the strait. Our study is the first to assess ship strike risk in this region of high whale density and marine traffic use.


INTRODUCTION
Vessel strikes, or collisions between ships and cetaceans, are a key threat to the recovery of baleen whale populations in many areas of the world, including Canadian Pacific waters (Fisheries and Oceans Canada 2013a,b).Baleen whales are at greater risk of being struck by ships than other marine mammals because of their large body size (Silber et al. 2010, McKenna et al. 2015).These species also spend extended periods of time at or near the sur-face, either feeding (Kot et al. 2014, Constantine et al. 2015) or recovering from the energetic demands of lunge-feeding at depth (Acevedo-Gutiérrez et al. 2002, Goldbogen et al. 2006), which increases their vulnerability (Laist et al. 2001).Furthermore, most baleen whales exhibit a limited ability to manoeuvre away from close-approaching vessels, or do not attempt to avoid ships at all (Nowacek et al. 2004, Harris et al. 2012, McKenna et al. 2015).This lack of avoidance behaviour may be caused by habituation to vessel noise (Nowacek et al. 2004), failure to per-ceive the vessel as a threat, or unwillingness to cease important activities such as feeding (Panigada et al. 2006, Silber et al. 2010).
Humpback Megaptera novaeangliae and fin whale Balaenoptera physalus populations in British Columbia (BC) are listed as 'Special concern' and 'Threatened', respectively, under Canada's Species at Risk Act (COSEWIC 2015).Both populations have been undergoing recovery after severe depletion by commercial whaling, which ended in the Canadian Pacific in 1967 (Ford 2014).Humpback abundance in coastal BC was estimated at 2145 individuals (ind.) in 2006, with an annual growth rate of about 4% (Ford et al. 2009).Qualitative impressions during field studies suggest that fin whale abundance in BC is also increasing, although likely not as rapidly as it has for humpback whales (Ford 2014).Vessel strikes have been identified as an important conservation concern that threatens the continuing recovery of both fin and humpback whale populations along the BC coast (Fisheries and Oceans Canada 2013a,b).During 2004−2011, 1 fin whale and 20 humpback whales were reported struck by vessels in BC (Fisheries and Oceans Canada unpubl.data).Humpback whales were the most commonly reported species involved in vessel collisions, with an individual reported injured or killed approximately every 9 mo.They were also the most frequently observed species bearing healed or partially healed wounds indicative of vessel collision injuries (Fisheries and Oceans Canada unpubl.data).
Due to the difficulty of recovering carcasses for necropsy and obtaining eye-witness reports, documented strike rates significantly underestimate the true impact of vessel collisions on whale populations (Ford et al. 2010, Conn & Silber 2013, McKenna et al. 2015).In particular, these sources of information about vessel strikes are biased toward near-shore areas and smaller vessels (versus larger, ocean-going cargo ships or tankers).Thus, strike rates for species with primarily offshore distributions, such as blue whales and fin whales (Ford 2014), go largely undocumented.Fin whales photo-identified in BC seldom bear scars attributable to ship strikes (Fisheries and Oceans Canada unpubl.data), which suggests that most individuals do not survive being struck (or that few fin whale strikes involve smaller vessels, which are more likely to cause non-fatal wounds) (Panigada et al. 2006).Strikes by very large ships are often not detected by mariners because collision impacts are unlikely to be felt and the bows of large ships are generally not visible to their crews (Conn & Silber 2013).Spatial models of ship strike risk based on the overlap between whale populations and marine traffic can complement information provided by necropsies and eye-witness reports.These models are able to predict ship strike risk over large areas and identify the regions of highest conservation concern where injuries are either fatal or compromise vital life processes, such as feeding and reproduction, with potential population-level impacts.Spatial analysis of variation in vessel speeds is critical to identifying areas where strikes are most likely to be lethal.Ship speeds exceeding 10 knots are more likely than not to cause mortality, and speeds ≥18 knots are almost certain to be lethal (Vanderlaan & Taggart 2007, Conn & Silber 2013).Information about strike risk offshore of the west coast of Vancouver Island is particularly vital because this area is a high-use region for marine traffic, especially for large, fast-moving commercial ships transiting Juan de Fuca Strait, a major shipping channel that provides access to several large ports (Vancouver, BC, and Seattle and Tacoma, WA).Between 10 000 and 11 000 vessels of all types enter this confined waterway every year (Nuka Research and Planning Group 2013).
Here, we perform the first spatially explicit analysis estimating the relative risk of lethal collisions between ships and whales off the west coast of Vancouver Island, Canada.We calculate strike risk for 2 of the most frequently observed species of baleen whales in BC, humpback and fin whales.We predicted whale densities across the study region using systematic aerial surveys (2012−2015), and overlaid these densities with a marine traffic dataset (2013) obtained from automatic ship tracking information collected by the Automatic Identification System (AIS) to calculate the relative risk of both vessel strikes and collision lethality.

Whale sightings data
Cetacean surveys were conducted from a De Havilland DHC-8-102 Dash-8 aircraft flown along systematically placed transects at a nominal speed of 278 km h −1 (150 knots) and an altitude of 305 m (1000 ft).Transects ran northeast to southwest, roughly perpendicular to the west coast of Vancouver Island at intervals of approximately 16 km.Two observers positioned at special large observation windows aft of the cockpit (left and right) reported all whale sightings to a data recorder, who entered them into a laptop computer.Sightings were reported as the whale(s) passed perpendicular to the aircraft, when observers measured an angle of declination to each sighting using a hand-held clinometer.Declination angles were reported to the nearest degree, such that 0° represented a sighting on the horizon, and 90°r epresented a sighting directly below the airplane.Angles between 70−90° were infrequently reported because this section of the water was generally not visible to the observers (who were scanning through flat, not bubble, windows) when the aircraft was flying on the level.Measurement error in reported declination angles (θ) affects the calculation of sighting distance from the transect (d).This is particularly true for distant sightings with small values of θ.However, measurement errors in θ likely did not exceed 5°, meaning that imprecision in distant sightings could have been on the order of several 100s of metres, which is not enough to affect the overall distribution of sightings, given the 25 km 2 grid cell resolution used in the subsequent models.
Once a sighting was reported, and if observers required additional time to identify species and number of individuals, the plane was flown in a loop around the whale(s).Once this was accomplished, the aircraft rejoined the survey transect.Observers also reported environmental conditions (sea state, visibility, precipitation, glare) using standardized categories at 5 min intervals throughout each survey, whenever conditions changed, and at the beginning and end of every transect.Geographic positions along the survey route were recorded automatically using the aircraft's GPS, at a sampling rate of either 1 or 0.2 Hz, depending on the survey year.
Effort and sightings data were filtered based on the recorded environmental conditions and survey status.Only 'on effort' sightings and survey track lines were included, and re-sightings were discounted from the analysis.Additionally, any 'on effort, closing' track lines, such as loops made by the aircraft to assist in species identification or group size counts, were excluded from the final effort data.The effort tracks and associated sightings that occurred during sea states > 4 (Beaufort wind force scale) or when visibility was reduced to ≤5 nmi (9.25 km) from the aircraft were also excluded.Occasionally, if observers could not positively identify a whale to the species level, but deemed it highly likely to be a particular species based on its morphology or behaviour, it was categorized as 'like humpback whale' or 'like fin whale'.These probable sightings were incorporated into the final sightings tallies used to model whale densities.
To estimate the true geographic position of each whale sighting, we used the following procedure.First, we calculated the perpendicular distance of every sighting from the aircraft using a formula from Buckland et al. (2001): where d is the distance (m) of the whale(s) from the transect, a is the altitude (m) of the aircraft, and θ is the declination angle (rad) formed between the horizon and the whale(s).We discounted sightings without reported declination angles or where θ = 0, because no horizontal distance could be calculated in these cases.We then found the compass bearing to each sighting by adding or subtracting 1.57 rad (90°) from the heading of the airplane (0 rad = north, and heading increased in a clockwise direction) at the time the sighting occurred, depending on whether the sighting was on the right or left side, respectively.Negative bearings and those > 6.28 rad (360°) were corrected by adding or subtracting 6.28 rad to obtain the equivalent angle.We then estimated the geographic position of each sighting based on its distance and bearing, as well as the plane's location: where ψ 1 and λ 1 are the latitude and longitude (rad), respectively, of the aircraft at the time the sighting was reported, ψ 2 and λ 2 are the latitude and longitude (rad) of the whale(s), d is the distance (m) from the transect to the whale(s), R is the radius of the Earth (6 371 000 m), and θ is the compass bearing (rad) to the whale(s).

Aerial survey effort
To determine the cumulative area (km 2 ) surveyed (and thus account for differences in the spatial distribution of effort), we calculated the width of the surveyed area for each transect and then summarized the variation in effort across a gridded surface of the study region.We began by constructing an effort buffer on both sides of every transect to determine the area that was effectively surveyed for whales.We excluded the section of the water directly beneath the plane, as it was not visible to the observers through the survey aircraft's flat windows.Given a reported average maximum sighting angle of θ = 1.22 rad (70°) below the horizon before the downward view became obstructed, and a nominal aircraft altitude of 305 m, we calculated this theoretical blindspot as follows (Buckland et al. 2001): and found it to have a total width of 222.6 m (i.e.111.3 m on either side of the transect line).Our effort buffer therefore excluded the strip extending from the transect line directly beneath the aircraft to a distance of 110 m on either side.We validated this theoretical blind strip by examining a histogram of the reported sighting distances, and found that sightings became extremely infrequent at distances <110 m from the transect line.
We determined the farthest extent of the effort buffer by constructing a detection function from the filtered, 'on effort' sightings of large baleen whales with the R package 'Distance' (Miller 2014) and calculated the resulting effective strip (half-) width (esw) (Buckland et al. 2001).We fit the preliminary detection function using conventional distance sampling (CDS) methods (Buckland et al. 2001, Thomas et al. 2002), with perpendicular distance as the only covariate.Prior to fitting the CDS detection function, 4 distance outliers (horizontal detection distances > 6000 m from the transect) were identified using Cleveland dot plots (Zuur et al. 2010) and removed from the dataset.Candidate detection functions included the hazard-rate and half-normal models, which were evaluated using Akaike's information criterion (AIC).Simple polynomial and cosine expansion terms were also considered.Left truncation was set at 1.25% in the initial detection function to ensure that all sightings <110 m from the transect (blindspot directly beneath the aircraft) were excluded.To test for possible effects of other covariates on whale detectability, we also applied multiple covariate distance sampling (MCDS).Since these additional covariates were either not significant (i.e. did not improve the model fit: 'cluster size', 'sea state', and 'visibility' covariates) or could not be included due to sample size limitations ('observer ID' covariate), we selected the CDS model with the lowest AIC value as the best-fit detection function.The right-truncation distance (w) was equivalent to the distance at which detection probability dropped below ~0.10, as recommended by Buckland et al. (2001).All sightings made at distances exceeding the truncation distance w were discounted from further analysis.We estimated the detection function goodness-of-fit by examining quantile-quantile plots and performing chi-squared and Kolmogorov-Smirnov tests.
We calculated effective strip (half-) width (esw, or μ) according to the following formula (Thomas et al. 2002): where P a is the probability that a randomly chosen animal within the surveyed area is detected, and w is the right-hand truncation distance of the detection function.We constructed the effort buffer such that its farthest extent was equivalent to the effective strip width (μ) (Gowan & Ortega-Ortiz 2014), since as many whales are detected beyond this distance as are missed within it (Thomas et al. 2002).
We built an effort buffer in ArcGIS (ESRI 2013) using the Geospatial Modelling Environment (GME) (Beyer 2012) that extended from the left truncation distance (110 m) to the esw (1010 m) on either side of the surveyed transects.We then divided the survey region into a grid of 25 km 2 cells and calculated the aggregate area surveyed per cell by summing the total area of overlapping effort buffers contained within each cell.Only grid cells containing survey effort (i.e.buffer area > 0 km 2 ) were retained for subsequent analysis (n = 1636).The whale sightings, corrected for geographic position and weighted by cluster (group) size, were then summed within each of these grid cells.

Vessel traffic data
We analysed the spatial distribution of marine traffic using AIS data collected by the Canadian Coast Guard (CCG) in 2013.AIS-equipped vessels broadcast information about their position, course, and speed over ground (SOG) using very high frequency (VHF) radio signals, at sampling rates of several times per minute.The 2013 AIS dataset consisted of vessels that were legally obligated to participate in the AIS network, as well as those that were voluntarily equipped with AIS.Compulsory reporting applied to all ships (other than fishing vessels) belonging to the following categories: ≥500 gross tons (GT), ≥300 GT that were transiting international boundaries, and ≥150 GT that were travelling internationally and carrying >12 passengers (Simard et al. 2014).Voluntarily-equipped vessels included fishing boats and AIS fishing beacons.Simard et al. (2014) compiled these AIS data and binned the resulting traffic densities (measured in daily ship-hours [ship-h], averaged over the entire year) into 5 categories of ship speed (2−5, 5−10, 10−15, 15−20, and > 20 knots) across a grid of 1 km 2 cells.Ship speeds were determined from AIS positions using a multi-step filter that excluded speeds > 40 knots and smoothed sudden changes in speed using a 900 s moving average (Simard et al. 2014).Vessels not underway and stationary AIS fishing beacons (i.e.SOG ≤1 knot) were discounted.Additionally, our analysis excluded the slowest traffic category (2−5 knots) reported by Simard et al. (2014), as vessels travelling at such low speeds are unlikely to pose a lethal strike risk to whales (Vanderlaan & Taggart 2007, Conn & Silber 2013).Marine traffic included in the 2013 AIS dataset can be generally categorized into the following types: cargo (e.g.container ships, bulk carriers), tanker, passenger (e.g.cruise ships, ferries), tug and towing, fishing, and pleasure vessels (Simard et al. 2014).The first 3 categories (cargo, tanker, and passenger) are of most concern when assessing lethal ship strike risk to whales, given the typically greater sizes and speeds of these vessels.More detailed information about AIS data collection and processing is provided by Simard et al. (2014).

Relative probability of a whale−vessel encounter
Determining the relative probability of a ship strike (using the proxy of a whale and a vessel occupying the same grid cell) requires estimates of the relative probability of encountering whales and encountering vessels across all grid cells (Vanderlaan et al. 2008).To accomplish this, we first estimated humpback and fin whale densities (ind.per 25 km 2 cell) from the aerial survey sightings using generalized additive models (GAMs).Modelling of whale densities was limited to these 2 species, as other baleen whale species were either not observed during the aerial surveys (e.g.North Pacific right whales and sei whales) or were sighted so infrequently that construction of a spatial model was impossible (e.g.blue whales, grey whales, and minke whales).Prior to model construction, we undertook data exploration following the protocol described by Zuur et al. (2010) to ensure that underlying model assumptions were not violated.One outlier was removed from the humpback count data because this grid cell contained a single sighting but had a very small surveyed area (0.03 km 2 ), resulting in a misleadingly low predicted density that substantially influenced the dispersion of the data set.Potentially nonlinear relationships between explanatory variables and whale counts were assessed prior to inclusion in the candidate GAMs by building separate generalized linear models (GLMs) and fitting GAMs to the GLM residuals for each covariate in turn (Zuur 2012).Variables that displayed nonlinear relationships with the GLM residuals (effective degrees of freedom, or edf > 1) were included as smoothers in the final negative binomial GAMs, whereas those with edf = 1 were included as beta terms.We constructed a set of candidate GAMs with latitude (converted to Universal Transverse Mercator [UTM] northing), slope, and depth as possible explanatory environmental variables (see Table 1).Longitude was excluded as an explanatory variable, as it was highly correlated with depth (r = 0.76, p < 0.0001) in our study region.An offset term to account for relative survey effort per cell (aggregate buffer area in km 2 ) was also included in all candidate models.
All survey effort and sightings were aggregated and analysed as a single dataset for each species.Both humpback and fin whales are present yearround in Canadian Pacific waters (Mizroch et al. 2009, Ford 2014); however, greater numbers of humpbacks occur from spring through fall (April− November) than during the remainder of the year.In contrast, there is little evidence of a distinct seasonal pattern to fin whale occurrence (Ford 2014).However, to ensure that aggregating the survey data across months did not impact model inference, we produced plots of GAM residuals (models run on pooled data) partitioned by month (Zuur 2012) for each species.The distribution of these residuals was similar across months, indicating that a factor covariate of 'month' would add little explanatory power to the models.
We fit negative binomial GAMs (logarithmic link function) in R using the 'mgcv' package (Wood 2004(Wood , 2011) ) at a spatial resolution of 25 km 2 for both species, and then used these models to predict whale densities across the study region at a finer resolution of 1 km 2 , to match that of the AIS ship data.Given that relationships between cetaceans and their habitats are scale-dependent, spatial scale is an important consideration when developing cetacean− habitat models (Redfern et al. 2006).Worldwide, satellite-tagging studies of fin and humpback whales have indicated that individuals travel at average rates of ~2−8 km h −1 , and even during area-restricted search (ARS) behaviour (presumed foraging), mean speeds typically exceeded 1.5 km h −1 (Silva et al. 2013, Kennedy et al. 2014, Rosenbaum et al. 2014).Changes in the environmental covariates (latitude, slope, and depth) at a 1 km 2 resolution are likely too fine-scale to noticeably influence the distribution of large rorquals that transit these distances easily in a fraction of an hour.We therefore chose to fit the GAMs to a larger resolution grid (25 km 2 ) than our ship data (1 km 2 ) because this scale was deemed more biologically relevant for predicting the distribution of whales.During model fitting, the appropriate smoothness for each covariate was estimated using likelihood-based methods (restricted maximum likelihood [REML]).We chose the negative binomial error distribution for fitting the GAMs, as the response variable (N ) consisted of over-dispersed (zero-inflated) count data, and global GAMs took the general form: (6) The depth covariate was square-root-transformed to make this predictor variable more uniform and thus reduce differences in the leverages of individual data points, which helped to stabilize model predictions (Wood 2006).We selected the best-fit GAMs for predicting humpback and fin whale densities by comparing the AIC scores (Zuur et al. 2009) of the various candidate models, which were generated from each global model by using backwards selection to drop non-significant covariates.GAM overfitting was avoided by incorporating the multiplier gamma = 1.4 (Kim & Gu 2004) to inflate the edf in the REML score.To assess whether or not spatial autocorrelation was present in the model residuals, we plotted variograms using the 'gstat' package in R (Pebesma 2004), and also examined the spatial distribution of residuals by size throughout study area.
We used GAM-predicted whale densities (W i ) from the top-ranked models to estimate the probability of observing a humpback or fin whale (P rel (Whale) i ) within each grid cell i, relative to all other grid cells (n), with an approach adapted from Vanderlaan et al. ( 2008): (7) P rel (Whale) i was calculated separately for fin whales and humpback whales.We likewise standardized the vessel traffic intensity values to determine the relative probability of observing a ship (P rel (Vessel) i ) within each grid cell i, over the total study area of n cells (Vanderlaan et al. 2008): (8) where V i was the annual average of daily ship-h km −2 (vessels travelling at speeds ≥5 knots) within each cell i.
The relative probability that a whale and a vessel encounter one another (i.e.occupy the same cell) was used as a proxy for the risk of a vessel striking a whale within each grid cell i in a domain of n cells (Vanderlaan et al. 2008): (9) where estimates of P rel (Encounter) i were also standardized such that their sum in the domain of n grid cells was equal to one.P rel (Encounter) i was calculated for fin whales and humpback whales separately.

Risk of a lethal strike as a function of ship speed
We determined the mean vessel speed (knots) per surveyed grid cell (Speed i ) using the following formula: (10) where was the median speed of each of the 4 vessel speed classes (7.5, 12.5, 17.5, and 23 knots, respectively), and Density c was the vessel traffic intensity (mean annual daily ship-h) per cell i for each speed class c.We then determined the probability that a ship strike would inflict a lethal injury using the mean vessel speed per cell and a logistic regression model developed by Conn & Silber (2013): (11) We then estimated the relative risk of a lethal collision between a vessel and a whale as follows ( Vanderlaan & Taggart 2007): To identify areas with the highest relative risk of collisions and collision mortality, we extracted the grid cells representing the 95 th percentile of P rel (Encounter) i and RR i values for both humpback and fin whales, and compared these regions to the remainder of the study area.Probabilities of whale, vessel, and whale−vessel encounters (i.e.ship strikes), as well as mean ship speed and relative risk of lethal whale−vessel encounters, were mapped for humpback and fin whales using the 'PBSmapping' package (Schnute et al. 2015) in R. All summary statistics are presented as mean ± SD, unless otherwise specified.

N s UTM northing s depth slope offset effo
) ) .

Whale sightings data and aerial survey effort
We conducted aerial surveys on 34 days from 2012−2015, during all months of the year except April, May, and August.Excluding those sections with poor environmental conditions, or where the plane flew in a closing loop or at altitudes > 366 m, we surveyed a total of 21 801 km of 'on effort' line transects.Since the majority of surveys took place in the fall and winter months, with the greatest number occurring in September (n = 11), our estimates of relative strike risk are therefore most applicable across these seasons.
The detection function with the lowest AIC value (4054.86)was a CDS hazard rate model with no expansion terms and a right-truncation distance (w) of 2650 m (Fig. 1).From this detection function, we calculated the effective strip (half-) width (esw, μ), or the farthest extent of the effort buffer, as 1010 m on either side of the aircraft.The aggregated effort buffers from all 34 surveys comprised a total surveyed area of 39 120 km 2 .After filtering for weather conditions, effort status, and altitude, 276 of the 322 total baleen whale sightings with distances less than the truncation distance (2650 m) remained.This included a total of 159 humpback whale or 'like humpback whale' sightings (329 ind.), and 74 fin whale or 'like fin whale' sightings (120 ind.; Fig. 2), which were input into the GAMs following the exclusion of a single outlier in the humpback sighting data.Mean group size per sighting was 2.1 ± 3.5 ind.for humpback whales (range = 1−33) and 1.6 ± 1.0 ind.for fin whales (range = 1−5; Fig. 2).There were also 3 sightings of single blue whales, all of which were observed west of the continental shelf break (Fig. 2).

Vessel traffic data
The 2013 AIS vessel data indicated that shipping traffic was less dense offshore, but became much more concentrated as it funnelled into or out of Juan de Fuca Strait (Fig. 3a).In particular, a commonly transited route is apparent that begins offshore (around 48.5°N, 128.0°W) and becomes more heavily used by ships as it moves eastward, toward the entrance of Juan de Fuca Strait and its Traffic Separation Scheme (TSS) lanes.Mean shipping intensity per grid cell for the entire study area was 0.006 ± 0.018 daily ship-h km −2 (range = 0−0.44),and the mean relative probability of encountering a vessel per cell, P rel (Vessel) i , was 4.2 × 10 −5 ± 12.3 × 10 −5 (range = 0−0.003).AIS-reporting traffic in 2013 (excluding the slowest 2−5 knot speed category and cells without traffic) travelled at mean speeds exceeding 12 knots throughout most of the study area = 12.5 ± 1.8 knots, range = 7.5−21.4knots; Fig. 3b).Vessel speeds were highest (≥16 knots) near the continental shelf break (200 m bathymetric contour) at the northern end of the study area, offshore of the shelf break, and inside Juan de Fuca Strait (Fig. 3b).Areas with average vessel speeds ≤10 knots were limited, and primarily occurred closer to Vancouver Island and at one location along the southern portion of the continental shelf break (Fig. 3b).Slow speeds at the shelf break location were likely the result of speed contributions from vessels engaged in fishing activities (see Simard et al. 2014, their Fig.13).

Relative probability of a whale−vessel encounter
The top-ranked GAM for predicting humpback whale densities included both latitude (UTM northing) and depth as explanatory variables, while the top-ranked GAM for fin whale densities only included depth (Table 1).Model averaging was not applied to the 2 highest ranked humpback-whale GAMs; we retained the simpler model with only the depth and latitude smoothers because the slope term was non-significant in the global model (p = 0.10).Model fits were significantly improved when depth was square-root-transformed. Detailed summaries of GAM outputs are presented in Tables 2 & 3  ranked GAMs predicted mean densities of 0.008 ± 0.014 ind.km −2 (max.= 0.089) for humpback whales and 0.003 ± 0.002 ind.km −2 (max.= 0.005) for fin whales.Effective degrees of freedom of smooth terms selected in the top-ranked models indicated that relationships between the significant environmental covariates and whale density were nonlinear (edf > 1; Fig. 4).Model smoothers for humpback whales predicted higher densities of individuals at the lowest latitudes in the study area (~48.1°N, or level with Cape Flattery and the Washington State coast) and at intermediate latitudes (~49.3°N, or level with Nootka Island and Hesquiat Peninsula on the Vancouver Island coast) (Fig. 4).The highest den-sities of humpback whales were predicted to occur in areas with depths of ~200 m (Fig. 4), which represents the edge of the continental shelf in our study region (the continental slope begins around the 200 m bathymetric contour and eventually levels out again at ~2300 m, where it reaches the abyssal plain).
The highest densities of fin whales were predicted to occur off the edge of the continental shelf, in water depths > 450 m (Fig. 4).Regions with the highest relative probabilities of encountering whales (per 1 km 2 cell) reflected the predictions of the GAM smoothers: humpbacks were most likely to be found along the continental shelf edge (200 m isobath) and in Juan de Fuca Strait (Fig. 5a), while fin whales were most likely to be found offshore of the shelf break (> 450 m; Fig. 5b).Fin whale encounter probability was very low in shore of the shelf break (Fig. 5b).Predictive power of the top-ranked humpback density model was fairly high, with 66.1% deviance explained, while the top-ranked fin whale model had a lower ex plained deviance of 26.7%.This is likely because only a single explanatory covariate (depth) was retained as significant in the top-ranked fin whale model.
The mean relative probability of a vessel encountering a whale (proxy for a ship strike) off the west coast of Vancouver Island was 4.2 × 10 −5 ± 23.8 × 10 −5 km −2 (max.= 0.007 km −2 ) for humpback whales and 4.2 × 10 −5 ± 6.3 × 10 −5 km −2 (max.= 0.001 km −2 ) for fin whales.Humpback whales were most likely to be struck along the continental shelf break, the   inshore approaches to Juan de Fuca Strait (east of 200 m isobath, ~48.5°N), and within the strait itself (Fig. 5c).The mean collision probability with humpback whales was 32.3-fold higher in these areas (95 th percentile) compared to the rest of the study domain.
Fin whales were most likely to be struck in the offshore approaches to Juan de Fuca Strait (west of the 200 m isobath, ~48.5°N) and inside the strait (Fig. 5d), where collision probability was 7.7-fold higher (95 th percentile) than the rest of the study domain.

Risk of a lethal strike as a function of ship speed
Across the study region, the mean relative risk of lethal ship strikes was 0.3 × 10 −4 ± 1.8 × 10 −4 km −2 for humpback whales (Fig. 5e) and 0.3 × 10 −4 ± 0.5 × 10 −4 km −2 for fin whales (Fig. 5f).In areas with the greatest risk of lethal ship strikes (95 th percentile; Fig. 6), we estimated the mean relative risk of lethal collisions with humpback whales to be 3.7 × 10 −4 ± 6.9 × 10 −4 km −2 , a 35.2-fold increase compared to the remainder of the study domain.These regions of highest relative lethal strike risk for humpbacks included Juan de Fuca Strait, an area due west of its entrance and inshore of the continental shelf break, and some areas overlying the 200 m bathymetric contour along the shelf break itself (Fig. 6a).For fin whales, the areas of highest concern had a mean relative lethal strike risk of 1.8 × 10 −4 ± 1.1 × 10 −4 km −2 (Fig. 6b), an 8.1-fold increase over the rest of the study region.Locations with the highest relative risk of lethal collisions with fin whales included Juan de Fuca Strait, as well as an area due west of its entrance (48.5°N) but offshore of the continental shelf break (Fig. 6b).Within the areas of highest lethal strike risk (Fig. 6), we estimated a mean value of the probability that a ship strike would be lethal due to speed (P(Lethal) i ) of 0.70 ± 0.09 for humpback whales and 0.75 ± 0.04 for fin whales, compared to a mean of 0.68 ± 0.10 (both species) in the remainder of the study area.In other words, a whale that is struck in these high-risk areas has a 70% (humpbacks) or 75% (fin whales) chance, on average, of being killed, compared to a 68% chance elsewhere in the study domain.

DISCUSSION
GAMs proved to be an effective and powerful method for estimating continuous, quantitative gra-dients of whale density across our study area using discrete point observations of individual whales or groups of whales.Spatial distributions of whale encounter probabilities estimated from these GAMpredicted densities corroborate existing information about humpback and fin whale distributions in BC.We determined that humpback whale distribution off southwestern Vancouver Island could be predicted based on latitude (UTM northing) and water depth, a Note that colour bar increments for the lowest and/or highest categories are not necessarily placed at equivalent intervals, to allow for improved visualization of the majority of the data range conclusion that is supported by similar models (Best & Halpin 2009, Williams & O'Hara 2009, Dalla Rosa et al. 2012) applied to humpback sightings over larger areas of the BC coast.Dalla Rosa et al. (2012) determined that humpbacks primarily favoured midshelf waters (50− 200 m), which matches our finding that humpback densities increased around the continental shelf break.Ford (2014) likewise describes humpback feeding areas in BC as primarily occurring in coastal or shelf waters, and a model based on historic whaling catch locations also indicated that humpbacks were predominantly distributed inshore of the shelf break (Gregr & Trites 2001).Although Dalla Rosa et al. (2012) predicted maximum humpback densities at slightly shallower depths (100 m) than our model (200 m), this is likely due to their larger study area that covered the entire BC coast.The shelf break in our smaller study domain off Vancouver Island occurs much farther from shore than in the rest of BC and the slope between the 100−200 m isobaths is also more gradual.This means that the 100 m isobath is not located near the shelf break off Vancouver Island, as it is in other coastal regions of BC (e.g.Haida Gwaii), and thus maximum humpback densities in our study area can reasonably be expected to occur in somewhat deeper water.In California, Dransfield et al. (2014) similarly predicted that hump backs had an affinity for the shelf break and occurred at greatest frequencies near the 200 m isobath.
As with our model, Dalla Rosa et al. (2012) determined that latitude was a significant predictor of humpback densities, and suggested this might occur because the extent of preferred, on-shelf habitat varies greatly by latitude in BC.Off Vancouver Island, the continental shelf ranges between 5 and 75 km wide, depending on latitude (Barrie et al. 2014).Higher humpback whale encounter rates are expected in regions where the continental shelf is wider, as more of the primary productivity is retained on-shelf in these areas, thus remaining available to coastal food webs (Perry et al. 1989, Ware & Thomson 2005).However, while the highest humpback densities in our study area were predicted at the lowest latitudes (~48° N) where the area of shelf habitat is also greatest, there was also an increase in estimated humpback density that occurred further north (~48.5°N), where the shelf is actually narrower.It is likely that gradients in densities of large whales with respect to latitude and/or depth are predominantly caused by underlying differences in prey availability in relation to these habitat variables (Dalla Rosa et al. 2012).
Our fin whale distribution model also matched that of Williams & O'Hara (2009) in that depth was a significant predictor of density; however, we did not find that fin whales were associated with latitude features.This difference is most likely due to the smaller size of our study area, which comprised a much narrower latitudinal range than that of Williams & O'Hara (2009).Fin whale sightings in BC have been primarily located in deep water beyond the continental shelf break (Ford 2014), which supports our model's predictions.Habitat models based on whaling catches for fin whales in BC also confirm that this species' historic distribution was largely offshore of the continental shelf (Gregr & Trites 2001).Although fin whales also regularly occur in certain inshore areas of the BC coast (Gregr & Trites 2001, Best & Halpin 2009, Ford 2014), none of these coastal habitats are located within our study domain.While quanti tative information about the specific distribution of fin whales by depth is not available for BC, fin whale distribution studies in the Mediterranean Sea also found that this species is primarily observed in offshore waters beyond the continental shelf (Forcada et al. 1996, Panigada et al. 2008).As in our model, the highest densities of fin whales in the Mediterranean oc curred at depths >1000 m, and declined sharply with decreasing depth (Panigada et al. 2008).
GAM-predicted mean densities of humpback whales (0.008 ind.km −2 ) were more than double that of fin whales (0.003 ind.km −2 ).This is consistent with recent assessments and surveys of cetacean abundance in BC, which report humpbacks as the most commonly observed species of baleen whale (Ford et al. 2010, Dalla Rosa et al. 2012, Ford 2014).These higher densities reflect the strong population recovery that humpbacks have undergone since the end of commercial whaling in BC (Ford 2014).The smaller relative population size of fin whales is reflected in the lower densities and probability of encounter for this species predicted by our model, as well as the lower frequency of fin whale sightings during the aerial surveys.We also found only a small degree of overlap between our model-predicted humpback and fin whale distributions, suggesting the possibility of habitat partitioning in regard to prey type and/or patch characteristics (Zerbini et al. 2006).Such partitioning is possible, since fin whales primarily consume euphausiid zooplankton, whereas humpback diet is more diverse and includes both zooplankton and small schooling fish (Flinn et al. 2002, Ford 2014).
Our model estimates of humpback and fin whale densities off Vancouver Island are conservative rela-tive to true whale densities, which are likely higher.This means that predictions of relative lethal strike risk calculated from the model-estimated whale densities should also be interpreted as conservative, minimum estimates.We were unable to correct for missed animals resulting from either availability bias (diving whales that were underwater, and thus unavailable for observers to detect), or from perception bias (whales that were available at the surface for observers to detect, but were not seen due to environmental conditions, fatigue, etc.) (Marsh & Sinclair 1989, Buckland et al. 2001, Dawson et al. 2008).For aerial surveys conducted in other regions, availability bias has been calculated for humpbacks at 0.67 (breeding ground off Brazil;Andriolo et al. 2006) and 0.26 (breeding ground off Hawaii; Mobley et al. 2001), and for fin whales at 0.25 (feeding ground in the Mediterranean Sea; Bauer et al. 2015).While we cannot apply measures of availability bias from other study areas to our calculations, these values give us a rough indication of the level of underestimation that may be present in our whale counts.Incorporating corrections for availability and perception bias would result in higher (and also more accurate) density estimates.Alternatively, it could be argued that diving whales, which are unavailable for counting, are similarly unavailable to be hit by a vessel.Actual densities of whales in danger of ship strike (i.e.near-surface) at any given time likely falls somewhere between the minimum whale densities we present here (uncorrected for biases), and the actual abundance of whales in the study area.Nonetheless, it is the spatial pattern and the relative differences in density between species (humpback and fin whales), and between various regions in the study area for each species, that provide the most important information for potential conservation actions.
Vessel traffic encounter probabilities in the study area increased toward the entrance of Juan de Fuca Strait, as this represents an area where ships enter or leave for the open sea and traffic approaching from the open sea becomes concentrated with increasing proximity to the designated shipping lanes (HEERC 2014).The route is regularly transited by large, deepsea cargo or tanker-type ships (HEERC 2014, Simard et al. 2014).These vessel types are of particular concern because of their large size (container ships are often > 300 m long with a beam of > 40 m; HEERC 2014), deep drafts (8−18 m; Silber et al. 2010, Constantine et al. 2015), and high speeds (typically 10−20 knots), which means that collisions with whales will most likely be lethal.Higher vessel speeds are correlated with shorter contact durations and increased accelerations experienced by struck whales, both of which lead to increased impact severity, and presumably, greater tissue damage (Silber et al. 2010).Mean vessel speeds were >12 knots throughout the majority of our study area off southwest Vancouver Island, which, in the event of a vessel strike, corresponds to a lethality probability of ~0.67 or higher (Conn & Silber 2013).In other words, throughout most of the study domain, whales struck by vessels would have > 67% chance of being killed.Higher vessel speeds are not only associated with greater strike mortality rates, but also lead to increased collision frequencies (Conn & Silber 2013, Lammers et al. 2013).In addition to direct strikes, large ships travelling at high speeds are also more likely to collide with whales as a result of hydrodynamic draw, which can pull a nearby whale toward the vessel's hull and thus extend the lethal strike zone to 1−2 times beyond a ship's actual draft (Silber et al. 2010).
We predicted the highest probabilities of whale− ship encounters (strikes) in regions where high whale densities co-occur with high-intensity maritime traffic, namely along the continental shelf break at the 200 m isobath (humpbacks), offshore of the shelf break (fin whales), and inside western Juan de Fuca Strait and the area west of its entrance (both species).Regions of highest risk for lethal collisions closely mirrored those with the highest incidence of strikes, since mean vessel speeds exceeded 12 knots throughout most of the study area.Even though fin whales occurred at much lower densities than humpback whales, the mean relative risk of a lethal collision per cell was actually quite similar for these 2 species (0.30 × 10 −4 km −2 for fin whales versus 0.29 × 10 −4 km −2 for humpbacks).This is likely related to the primarily offshore distribution of fin whales, which exposes them to marine traffic travelling at higher speeds, and thus results in equally high probabilities of lethal injury in the event of a strike, despite their lower densities.
We did not account for species-specific differences in vulnerability to ship strikes, and for fin whales especially, this may mean that risk is under-estimated.Worldwide, fin whales are the most frequently struck by vessels (Laist et al. 2001).Fin whale dive behaviour or anatomy may make them more vulnerable to lethal vessel strikes than other species, which our model did not take into account.Behavioural and physical factors that increase the vulnerability of some cetacean species to ship strikes include longer surface intervals, larger body size, limited ability to avoid approaching vessels due to either swim speed or manoeuvrability constraints, and auditory limitations such as directionality and distance perception (Lawson & Lesage in press).From a conservation standpoint, ship strike mortalities may also have a proportionally greater impact on the fin whale population because this species has a lower overall abundance than humpback whales, and is of greater conservation concern (listed as Threatened under the Species at Risk Act; COSE -WIC 2015).Estimates of potential biological removal (PBR) give an indication of the threshold above which human-caused mortality affects population recovery.PBR is an estimate of the maximum number of animals (excluding natural mortality) that may be removed per year, while still allowing the population to reach or sustain its 'optimum sustainable population' size (Wade 1998).PBR estimated for humpback whales in the Canadian Pacific, based on a population estimate in 2006 of 2145 ind.(range = 1970−2331), was 21 whales yr −1 (Ford et al. 2009).It is unlikely that a PBR estimate for fin whales would allow for more mortality than this, given the overall comparatively lower abundance of this species in BC (Ford 2014).
Despite having very low predicted encounter probabilities for fin whales, and only moderate encounter probabilities for humpbacks, the western portion of Juan de Fuca Strait is a relatively high-risk area for lethal collisions with both species.This is likely due to a combination of factors: the strait has a very high intensity of vessel traffic and higher than average vessel speeds.Therefore, any whales found within Juan de Fuca Strait, however infrequently, are exposed to high risk of a lethal collision.In addition, GAMs predict non-zero densities of whales across every cell in a study domain, not just those in which actual sightings were recorded -so cells where fin whales are unlikely to occur will actually show positive (albeit very small) density estimates, which translate into high whale−ship encounter probabilities in areas where vessel traffic is very intense.Although no fin whales were observed inside Juan de Fuca Strait during our aerial surveys, infrequent sightings of fin whales have been reported (Ford 2014, Chamberlain 2015, Mark Malleson pers. comm. November 2015), indicating that fin whales do occasionally enter the strait.For this reason, we believe the predicted increased strike risk to fin whales within the strait is reasonable, in the rare instances that they are found in this location.In other studies, physical bottlenecks where marine traffic becomes concentrated (like Juan de Fuca Strait) have been similarly associated with increased strike risk (Williams & O'Hara 2009, Silber et al. 2010).
Strike rates could be higher at times of the year when whales are most abundant in an area, or when marine traffic is highest, if indeed there is seasonality to commercial shipping (Panigada et al. 2006, Lammers et al. 2013).Our analysis was based on an aggregate dataset of all aerial survey sightings (2012−2015) and all AIS vessel traffic (2013), regardless of season, and our whale distribution models contained only 1 or 2 static explanatory variables.Future models could improve upon both the spatial and temporal resolution of predicted whale densities by incorporating dynamic, time-varying environmental predictors, as well as additional further survey data to support finer-resolution models.
Our analysis of the spatial distribution of ship strike risk off the west coast of Vancouver Island provides the first step towards ascertaining the potential need for mitigation strategies to reduce the threat of lethal vessel strikes to humpback and fin whales in this area.In other regions worldwide, where the risk of ship strikes has been found to be a significant threat to species conservation or recovery, the most commonly used and effective vessel strike mitigation strategies generally fall into 2 categories: speed limits for ships transiting areas with high whale densities to reduce both the rate and lethality of collisions (Wiley et al. 2011, McKenna et al. 2012, Conn & Silber 2013), or diverting traffic to avoid such areas entirely and thus reduce the co-occurrence of ships and whales (Vanderlaan & Taggart 2009).Traffic Separation Schemes (TSS) that intersect important whale habitat can also be repositioned to help prevent vessel strikes (Vanderlaan et al. 2008, Silber et al. 2012).Other conservation strategies that have been implemented to reduce ship strikes include passive acoustic monitoring (PAM) buoys or bottommounted arrays, which detect whale vocalizations in near real-time and broadcast alerts using vessel communication technologies such as AIS or Navigational Telex (NAVTEX) (Clark et al. 2007, Van Parijs et al. 2009, Morano et al. 2012, Reimer et al. 2016), as well as mandatory VHF marine radio reporting systems (Ward-Geiger et al. 2005, Brown et al. 2007, Silber et al. 2012).Both these approaches serve to warn mariners entering areas with high whale densities, but may be most effective with species that vocalize regularly, such as the North Atlantic right whale.Autonomous underwater vehicles (AUVs, or ocean gliders) have also been suggested as a tool for reporting real-time acoustic detections of whales to vessel operators (Baumgartner et al. 2013).Placing dedicated observers on vessels transiting high-risk areas might also reduce strike risk by increasing the likeli-hood that whales are detected in time for the ship's crew to take evasive action.Active acoustic alarms intended to warn whales of approaching vessels have been tested but were unsuccessful, as these sounds caused diving right whales to return to the surface, which increased their exposure to ship strikes rather than reducing it (Nowacek et al. 2004).
Our models provide the first spatially explicit information that can be applied to the development of conservation policies to reduce ship strike risk to humpback and fin whales off the west coast of Vancouver Island.In particular, discussions with the International Maritime Organization (IMO) and shipping industry could focus on the feasibility of reduced ship speeds or traffic re-routing in areas of highest concern: the continental shelf break off the west coast of Vancouver Island, and Juan de Fuca Strait.To improve upon the risk reduction strategies already suggested here, further surveys are needed to improve the temporal resolution of whale and ship distributions (i.e. to determine if there is a need for seasonally implemented mitigation strategies).In addition, the distributions of other cetacean species should be assessed to determine the potential impact that any mitigation efforts might have on these populations.

CONCLUSIONS
Ship strike risk analysis based on spatial models of whale density provided an effective approach for identifying regions of conservation concern, particularly in areas (e.g.offshore) where actual mortality rates are difficult to quantify from carcass evidence or eye-witness reports.We found evidence of habitat partitioning between humpback and fin whales, with humpback densities being highest on the continental shelf (particularly over the shelf break), and fin whales being distributed at somewhat lower densities in deeper, offshore waters.Whales had the potential to be struck by ships in any location where their distributions overlapped with that of marine traffic.In regions of both high whale densities and high-intensity marine traffic (e.g.Juan de Fuca Strait and the region due west of its entrance), whales were susceptible to elevated risk of lethal ship strikes.Ship speeds throughout the offshore area of the west coast of Vancouver Island were sufficiently high (>12 knots) that collisions with whales are more likely than not (> 50%) to result in lethal injuries.
The models we have developed for predicting the spatial distribution of vulnerable whale populations could be used to advise mariners about high-risk locations for ship−whale collisions.Ultimately, mitigation efforts to reduce the impact of ship strikes on whale populations could include speed restriction zones, areas to be avoided, or PAM-linked mariner notification systems (or a combination of these strategies).In addition, the modelled relative risk of lethal collisions could inform managers about the potential population-level impacts of ship strikes on humpback and fin whales off the west coast of Vancouver Island.Additional survey effort off the west coast of Vancouver Island could help to further refine our whale density models and ship strike risk estimates -for instance, more surveys might allow for analysis of seasonal or annual trends, or predictions for less-frequently encountered species (e.g.blue whales).

Fig. 3 .
Fig. 3. Colours indicate (a) probability of encountering a vessel and (b) estimated mean vessel speed (2013 Automatic Identification System [AIS] ship traffic dataset), in the aerial survey study area off the west coast of Vancouver Island, divided into 1 km 2 grid cells (n = 23 996).Filled cells indicate those areas containing survey effort that were retained for analysis.Continental shelf break indicated by the 200 m bathymetric contour (black line).Colour bar for (a) is scaled similarly to Fig. 6c,d to facilitate comparisons.Note that colour bar increments for the lowest and highest categories are not necessarily placed at equivalent intervals to allow for improved visualization of the majority of the data range

Fig. 4 .
Fig. 4. Smoothing functions (solid lines) with 95% CIs (shaded bands) for the explanatory variables, UTM northing (latitude) and depth, of the top-ranked negative binomial generalized additive models (GAMs) estimating (a,b) humpback Megaptera novaeangliae and (c) fin whale Balaenoptera physalus densities over a gridded surface of 25 km 2 cells.The y-axis labels display the fitted function with effective degrees of freedom in parentheses, while x-axis rug plots indicate the distribution of sampled values within each explanatory variable.For ease of interpretation, the latitude smoother includes x-axes showing both projected (UTM) and approximate geographic (latitude) values, and the depth smoothers include x-axes showing both square-root-transformed (m 1/2 ) and untransformed (m) values

Fig. 5 .
Fig. 5. Colours indicate relative probability of (a,b) encountering a humpback Megaptera novaeangliae or fin whale Balaenoptera physalus (calculated from generalized additive model [GAM] estimates of whale densities), (c,d) a vessel encountering a humpback or fin whale, and (e,f) relative risk of a lethal collision between a vessel and a humpback or fin whale (2013 Automatic Identification System [AIS] ship traffic dataset), in the aerial survey study area off the west coast of Vancouver Island, divided into 1 km 2 grid cells (n = 23 996).Filled cells indicate those containing survey effort that were retained for analysis.Continental shelf break is indicated by the 200 m bathymetric contour (black line).Colour bars for (c,d) and Fig. 3a are scaled similarly to facilitate comparisons, as are colour bars for (e,f).Note that colour bar increments for the lowest and/or highest categories are not necessarily placed at equivalent intervals, to allow for improved visualization of the majority of the data range

Fig. 6 .
Fig. 6.Filled cells (n = 1200) indicate areas of highest relative risk (95 th percentile) of lethal collisions (per 1 km 2 ) between ships and (a) humpback Megaptera novaeangliae and (b) fin whales Balaenoptera physalus; 2013 Automatic Identification System (AIS) ship traffic dataset.Mean relative risk of a lethal collision is 35.2-fold higher for humpbacks and 8.1-fold higher for fin whales in the illustrated areas than in the remainder of the surveyed study area.Continental shelf break is indicated by the 200 m bathymetric contour (black line).Colour bars are scaled similarly to one another, and to Fig. 5e,f to facilitate comparisons.Note that colour bar increments for the lowest and/or highest categories are not necessarily placed at equivalent intervals, to allow for improved visualization of the majority of the data range

Table 2 .
. The top- Negative binomial generalized additive model (GAM) (log-link function) summary of top-ranked model for humpback whale Megaptera novaeangliae counts (n = 159 sightings, 329 ind.) over the 25 km 2 gridded study area off the west coast of Vancouver Island.*Significant relationships (p < 0.05).AIC: Akaike's information criterion; edf: effective degrees of freedom AIC ΔAIC w Humpback whale candidate GAMs N ~ s (Ί |depth|) + a Final selected model

Table 1 .
Candidate generalized additive models (GAMs) for predicting densities of humpback Megaptera novaeangliae and fin whales Balaenoptera physalus off the west coast of VancouverIsland (2012Island ( −2015)).All models included the offset term of logged aggregate effort area (km 2 ) per grid cell.
w: model weight.Where candidate models had Akaike's information criterion (AIC) values differing by < 2, the more parsimonious model was chosen