Gray whale densities during a seismic survey off Sakhalin Island, Russia

Some whale populations that were severely reduced by commercial whaling have shown strong recovery since becoming protected, while others remain depleted and of high conservation concern. Small populations are particularly susceptible to anthropogenic threats, in cluding acoustic disturbance from industrial activities such as seismic surveys. Here, we investigated if sound exposure from a 16 d seismic survey displaced gray whales Eschrichtius robustus from their coastal feeding area off northeastern Sakhalin Island, Russia. We conducted multiple shore-based surveys per day, weather permitting, and created daily 1 km2 density surfaces that provided snapshots of gray whale distribution throughout the seismic activity. A Bayesian spatiotemporal analysis was used to examine possible effects of characteristics of sound exposure from seismic airguns on gray whale occupancy and abundance. Models suggested highest occupancy in areas with moderate sound exposure. Slightly decreased densities were associated with sound exposure when the pattern for the previous 3 d was high sound on Day 2 and low sound on Days 1 and 3. Our findings should be interpreted with caution, given the low number of positive densities. This was due to success of the primary mitigation measure, which was to conduct the seismic survey as early in the feeding season as possible when few gray whales would be present. It is also possible that observed differences in occupancy and density reflect changes in prey availability rather than noise. Prey distribution and abundance data were unavailable for our study, and this important covariate could not be included in models.


INTRODUCTION
Commercial whaling reduced many baleen whale populations to small proportions of their former abundance (Clapham et al. 1999).Protective measures under international law to stop commercial whaling of right whales Eubalaena glacialis and bowhead whales Balaena mysticetus were introduced in 1935, followed by gray whales Eschrichtius robustus in 1946, and humpback Megaptera novae angliae and blue whales Balaenoptera musculus in the mid-1960s (Best 1993).These measures culminated with the moratorium on commercial whaling in 1986.Some whale populations that were severely reduced by commercial whaling (e.g.see Clapham et al. 1999) have shown strong recovery since becoming protected (e.g.Calambokidis et al. 2008, Givens et al. 2010, Punt & Wade 2012).However, other populations, such as all North Atlantic and North Pacific right whale populations, and gray whales in the western North Pacific, remain depleted and continue to be of high conservation concern (e.g.Kraus et al. 2005, Weller et al. 2012).These small populations are particularly susceptible to anthropogenic threats such as entanglement in fishing gear, exposure to contaminants, habitat degradation, ship strikes, disturbance from vessel traffic, and hearing injury and acoustic disturbance from underwater sound (Clap ham et al. 1999, Reeves et al. 2003).
Seismic surveys that are conducted to map offshore oil and gas reserves produce substantial levels of low-frequency pulsed sound that can propagate over long distances (NRC 2005).Potential effects on marine mammals include permanent or temporary hearing injury, and behavioural disturbance (Nowacek et al. 2007, 2015, Southall et al. 2007).Mitigation to prevent hearing injury is now standard practice for seismic surveys in many parts of the world (Weir & Dolman 2007, Compton et al. 2008, JNCC 2010).However, disturbance is also of concern because recent work has postulated linkages between ob served behavioural changes in marine mammal species and population-level effects that are mediated through changes in life history functions and vital rates (NRC 2005, New et al. 2013, 2014).Indeed, preventing population-level effects on species of conservation concern is an important consideration when planning a seismic survey that can repeatedly en sonify an area they inhabit over several days.Southall et al. (2007) suggested that repeated or sustained disruption of life functions is more likely to affect vital rates than an isolated and brief disturbance, and that avoidance by the animals of important habitat can be significant if this effect lasts over multiple days or is recurrent.
Baleen whales are low-frequency hearing specialists (Southall et al. 2007), which means their hearing ranges overlap with the dominant frequencies produced by seismic airgun arrays.This makes baleen whales more sensitive to sound exposure and potential disturbance from seismic surveys compared to other marine mammal species (Nowacek et al. 2007).Numerous studies (reviewed in Richardson et al. 1995, Gordon et al. 2003, Nowacek et al. 2007) have reported acoustic disturbance in baleen whales that includes changes in respiration, movements and vocal behaviour, and avoidance of the sound source through short-or long-term displacement from ensonified areas, although considerable variability in species' responses has been noted.Recent literature suggests that not only total sound exposure levels, but also context (e.g.individual's age, sex, activity, habituation) influences the probability and type of behavioural response to anthropogenic sound (Wartzok et al. 2003, Southall et al. 2007, Ellison et al. 2012).Reported sound exposure levels associated with observed behavioural responses have ranged from approximately 120 to 180 dB re 1 µPa root mean square (rms) (reviewed in Gordon et al. 2003, Nowacek et al. 2007).For example, Malme et al. (1984Malme et al. ( , 1986) ) found that 10% of feeding gray whales avoided seismic airgun sound levels ≥ 163 dB re 1 µPa rms.
Sakhalin Energy Investment Company conducted a seismic survey from 17 June to 2 July 2010 on the northeast Sakhalin shelf, Russian Federation, adjacent to the southern portion of the primary known nearshore (Piltun) feeding area of gray whales in the western North Pacific.Acoustic propagation modelling of the 2620 in 3 (42 900 cm 3 ) airgun array used in the seismic survey predicted that sound levels greater than 163 dB re 1 µPa rms would occur in parts of the Piltun feeding area.A monitoring and mitigation plan was subsequently developed by the company in co-operation with IUCN's Western Gray Whale Advisory Panel to prevent behavioural disturbance to feeding gray whales in addition to standard miti gation for acoustic injury (Bröker et al. 2015).The primary mitigation measure for the seismic survey was to commence as soon as possible after ice melt; this timing coincided with the beginning of gray whale migration into the feeding area, when few animals would be present.
Standard statistical methods are often used to as sess if population-level avoidance of sound exposure from a seismic source occurred.These methods evaluate differences in mean counts or densities of a species or species group in selected areas when the airguns are active versus inactive.These areas may be within visual range of marine mammal observers on the seismic vessel (e.g.Stone & Tasker 2006, Weir 2008) or at distances farther from the seismic source using an independent survey platform (e.g.Richardson et al. 1999, McCauley et al. 2000).A more robust approach uses spatio-temporal analyses to explore changes in a species' spatial pattern over time and optionally to investigate covariate effects (e.g.Ver Hoef & Jansen 2007, Kirkman et al. 2013).While several methods are available for such analyses, regression using generalized linear mixed models (GLMM) within a Bayesian hierarchical modelling framework can be used to control for confounding covariates and address data depend-encies (Zhao et al. 2006).These dependencies include temporal autocorrelation from repeated sampling of spatial data (Bence 1995) and spatial autocorrelation in which data points closer in space have either more or less similar values than those farther apart (Legendre 1993).Species of conservation concern pose a particular challenge in spatio-temporal analyses because their low abundance results in many zero counts at sampling locations that cannot be modelled using standard statistical distributions (Cunningham & Lindenmayer 2005).Instead, zero-inflated models can be used that explicitly ac count for extra absences in data by modelling species abundance as 2 processes: (1) species presence, and (2) species abundance when present (Wenger & Freeman 2008).
The primary objective of the present study was to examine whether sound exposure from the 2010 seismic survey altered the distribution and abundance of gray whales within their feeding habitat.We conducted multiple shore-based gray whale surveys per day, weather permitting, and used these data to create fine-scale daily 1 km 2 density surfaces that provided snapshots of gray whale distribution in the southern part of the Piltun feeding area (see Fig. 1) throughout the seismic activity.These surfaces were used in a spatio-temporal analysis to assess the effects of total sound exposure levels and of characteristics of sound exposure from the seismic airguns on gray whale occupancy and abundance within grid cells of the density surfaces.A particular challenge in our study was the need to relate continuously recorded sound le vels from the seismic airguns to the daily (weather permitting) estimated gray whale grid cell densities.We employed novel methods, including functional principal components ana ly sis (FPCA; Ramsay & Silver man 2005), to convert time series of sound estimates into covariates that matched the temporal scale at which our res ponse of grid cell densities was estimated.

Study area and experimental design
The seismic vessel sailed 35 predefined lines spa ced ~300 m apart within an ~170 km 2 area (Fig. 1).Each line required ~2 h to complete at an average speed of 5 knots (9.3 km h -1 ).The vessel towed two 2620 in 3 (42 900 cm 3 ) airgun arrays that each consisted of 27 airguns configured in an alternate firing mode.Airguns were fired at ~8 s intervals.The airgun array was shut down during turns between lines, which took ~4 h.More details of the seismic survey are provided in Bröker et al. (2015).
Our study area was located in the southern portion of the Piltun feeding area that was closest to the seismic survey activity (Fig. 1).The ~50 km long by ~10 km wide study area extended from shore to a depth of shown as dotted blue and solid blue lines respectively ~20 m.Gray whales show high annual site fidelity to the Piltun feeding area, which has a high biomass of preferred benthic and epibenthic invertebrate prey (Fadeev et al. 2012) that varies spatially within the feeding area within and across years (Fadeev 2013).Although gray whale densities show considerable spatial and temporal variability in the feeding area, most whales are found within 5 km of shore, with peak densities between 500 and 2000 m from shore (~5 to 15 m water depths) (Vladimirov et al. 2013).Higher numbers of gray whales are typically found in the northern part of the study area that is adjacent to the lower half of the Piltun Lagoon (Vladimirov et al. 2013).
Our experimental unit was a survey of the study area that was used to estimate a 1 km 2 gray whale density surface (see next subsection for details).We conducted multiple shore-based surveys per day, weather permitting, to create a time series of density surfaces for analysis.As described in Muir et al. (2015a), 2 'distribution' teams jointly conducted scans from 5 onshore stations (9 to 13; Fig. 1) so that a complete survey of the study area was performed in approximately 2 h.These surveys were timed to coincide with either the sailing of a line when airguns were fired ('on' sound), or a turn between lines when the airguns were shut down ('off' sound).Distribution team effort was supplemented by 2 'behaviour' teams that conducted theodolite tracking and focal animal follows at the Seismic North and Seismic South stations in the central part of the study area (Fig. 1).Hourly scans were performed by a behaviour team throughout the day unless the team was engaged in focal-follow observations of a gray whale.
The scan protocols presented in detail in Muir et al. (2015a) are briefly described here.Observers used Fujinon FMTRC-SX 7 × 50 reticle binoculars to scan the nearshore waters surrounding a station at a constant rate (10°min −1 for distribution teams; 9.3°min −1 by behaviour teams).Bearing, reticle estimate, and number of individuals were recorded for each observed gray whale sighting during a scan.Environmental data (e.g.visibility, glare, Beaufort wind force scale [hereinafter Beaufort scale]) were collected for each scan at a station.Visibility was categorized as a 4 level code: (1) excellent conditions with clear horizon line, (2) good conditions with little to no haze and/or rain with relatively clear horizon line, (3) fair with some haze and/or rain but horizon still visible enough for reticle estimation, (4) poor, no visible horizon due to fog and/or rain.Scan surveys were not conducted if the visibility code was 4, Beaufort Scale exceeded 3, or wind speed exceeded 10.0 m s −1 .

Gray whale sighting data preparation, mapping and density calculations
Sighting locations were calculated and mapped as described in Muir et al. (2015a).The methodology to estimate density surfaces for shore-based scans is presented in Vladimirov et al. (2011) and summarized here.The study area was overlaid with a grid of 1 × 1 km cells, and a gray whale density surface was estimated for each survey.Effort and sightings for a survey included scans at all sampled distribution stations and any behavioural scans that overlapped the survey start and end time.We limited the radius of a scan's sampled area to the estimated 0.1 reticle distance (mean: 5245 m, range: 3907 to 6436 m) for the observer eye height at the scan's station.Sightings and effort beyond a station's 0.1 reticle distance were excluded.A density was estimated for each grid cell with at least 0.1 km 2 survey coverage.Some grid cells were sampled more than once during a survey by scans at adjacent stations (e.g.Stns 11 and Seismic South) or because multiple scans occurred at a behaviour station.Since scans were conducted at different times, they had slightly different patterns of 'on' and 'off' sound in the time preceding each scan.We therefore only allowed 1 scan to estimate a grid cell density per survey to ensure clear transitions in sound on/off patterns for grid cell sound covariate estimation (see 'Materials and methods: Sound co variates').Scans by distribution teams that worked throughout the study area were preferentially retained over behaviour scans.A single scan was randomly selected if adjacent distribution stations both sampled a grid cell.
A gray whale density (D ˆi,j ) was estimated in the j th grid cell that was sampled during survey i as follows: (1) where a ˆi,j is the estimated availability correction for the platform (distribution or behaviour) that sampled grid cell j during survey i (McLaren 1961 1 ); A i,j is the area covered by survey i in the j th cell; S i,j is the number of sightings by survey i in cell j; n i,j,k is the number of gray whales observed in the k th sighting by survey i in cell j; and p ˆis the probability of detecting an available gray whale.Previous work that analysed While the slightly modified method by Laake et al. (1997) may be more appropriate, use of McLaren (1961) here will not affect the results because the calculated availability correction was very similar for each of the behaviour and distribution densities (due to the slight difference in scanning rates) and was a common correction within each of these data sets gray whale sightings from a double-platform (vessel and shore-based) experiment estimated a flat detection function up to the 8 km distance tested (E.Rexstad & D. Borchers unpubl.),thus p ˆwas set to 1. Effects of environmental covariates on detection probability were not tes ted in the double-platform analysis due to small sample sizes.A density of zero was assigned to covered cells with no observed gray whales during a survey.Additional details of the density calculations, including methods to estimate the availability and detection probabilities, are provided in Supplement 1 (www.int-res.com/articles/suppl/n029p211_supp.pdf).

Sound covariates
Underwater sound generated by the airguns was measured throughout the seismic survey using 9 tele metered autonomous underwater acoustic recorders (T-AUARs) (Bröker et al. 2015, Racca et al. 2015).The T-AUARs were deployed along a 20 km section of the estimated Piltun feeding area boundary (Muir et al. 2015b) that was adjacent to the seismic survey area.These measurements of sound pulse levels were used in combination with numerical modelling to calculate a time series of 5 min cumulative sound exposure level (cSEL) bins through out the seismic activity at each grid cell centroid location.The cSEL calculations used the es timated median sound exposure level (SEL) over depth for individual seismic acquisition pulses and did not include contributions of non-pulse sound energy from other sources.Details of the acoustic numerical modelling are provided in Racca et al. (2015).
It is a frequent challenge in ecological studies to relate a continuously monitored physical measurement to a biological process that is infrequently sampled (Ainsworth et al. 2011).We used 2 approaches, summary statistics and summaries of sound patterns to convert the time series of airgun sound estimates at each grid cell to covariates that matched the temporal scale of our estimated survey-level grid cell densities.These covariates included total sound exposure and characteristics of sound exposure for each sampled grid cell in each survey's density surface over 3 time windows: short (3 h), intermediate (3 d) and long (7 d).A survey's grid cell time window ended at the midpoint time of the scan used to estimate that cell's density.
Seven types of sound covariates were calculated for each sampled grid cell using the cSEL time series in each time window for each scan (Table 1 The first 3 scores were used (3hPC1, 3hPC2, 3hPC3) FPCA 3 d scores FPCA scores calculated for 3 d time series using an 8 h smooth.

h smooth
The first 3 scores were used (3d8hPC1, 3d8hPC2, 3d8hPC3) (1 st 3 FPCs) FPCA 3 d scores FPCA scores calculated for 3 d time series using a 24 h smooth.

h smooth
The first 3 scores were used (3d24hPC1, 3d24hPC2, 3d24hPC3) (1 st 3 FPCs) Table 1.Grid cell sound covariates.Each time window used to calculate a covariate for a scan's grid cell ended at that scan's midpoint time.cSEL: cumulative sound exposure level, FPCA: functional principal components analysis fication level was determined by assigning each cell in the density surface grid to 1 of 3 categories (high, medium or low) based on the magnitude of the cell's 3 d cSEL cell value on Day 10 of the seismic survey (Fig. 2).The ensonification pattern covariate captured 3 different shapes of changing cSEL at a cell as the seismic vessel sailed by while firing the airgun array; these patterns depended on the cell's relative position to the seismic survey area and the vessel's sailing direction.We used FPCA (Ramsay & Silverman 2005) to summarize the patterns in cSEL across each of the 3 h, 3 d and 7 d time series at all grid cell centroids.
FPCA is analogous to standard principal components analysis (PCA) for reducing data dimensionality (Ramsay & Silverman 2005).While PCA reduces the number of variables, FPCA summarizes the pattern of data over time into a few variables (functional principal component scores [FPCs]).FPCA was run using the fda package in R v.2.15.2 (R Development Core Team 2012).The survey-level cSEL series in 3 h time windows were readily re presented as curves using bsplines.However, the 3 d and 7 d time windows spanned the sailing of several seismic lines, resulting in multiple jumps be tween on and off sound that were not well captured by bsplines.Instead, we smoothed the 5 min time series for the 3 d and 7 d time windows using 8 h averages, and again with 24 h averages, before calculating the FPCA covariates.The 2 sets of FPCs for each time window thus captured differences in the shape of the 8 h or 24 h means in sound within the time window.Details of the sound covariate calculations are provided in Supplement 2 (www.int-res.com/articles/ suppl/ n029p211_supp.pdf).

Model covariates and temporal blocking
Models included environmental co variates that could affect detection, and covariates for predicting space/ time changes in occupancy and abundance (Table 2).We used the Near tool in ArcGIS v.10.1 (ESRI 2012) to calculate perpendicular distance from shore to each grid cell centroid.Depth to the nearest metre was associated with each centroid using bathy metry generated from field sampled data (Caslys Consulting unpubl.).
After exploratory analyses of different time blocks (see Supplement 3 at www.int-res.com/articles/suppl/ n029p211_supp.pdf),we concluded that dai ly temporal blocking was needed to address challenges of (1) incomplete survey-level density surfaces due to poor weather or time constraints that occasionally prevented scans from being conducted at one or more stations during a survey, and (2) the high proportions of zero densities per surface.Mean daily densities were calculated in each sampled grid cell using the arithmetic average of that day's survey-level density estimates.Sound covariates within each grid cell were calculated for each daily density surface using the arithmetic mean of that day's survey-level cSEL covariates in dB (cSEL3h, cSEL3d, cSEL7d) and the FPCA scores.The proportion of non-zero sound, and median and standard deviation in 'on' sound, were calculated using all 5 min cSEL time series for the surveys being averaged.

Model structure and fitting
Zero-inflated Bayesian hierarchical GLMMs with conditional autoregressive random (CAR) effects (Ver Hoef & Jansen 2007) were used to model the effect of sound on gray whale densities.These models addressed several analytical concerns that included spatial correlation, a large proportion of zero densities, repeated sampling and increasing number of whales as they migrated into the study area.Zeroinflated models are decomposed into 2 parts: occupancy (presence/absence of gray whales in a grid cell) and abundance (gray whale density within occupied cells).A model can be formulated as either a '2-part' (or 'hurdle') model (Heilbron 1994), or a 'mixture' model (Lambert 1992).As its name implies, the 2-part model consists of 2 separate models: a logistic regression model of occupancy, and a zerotruncated regression model of abundance given presence.The mixture model differs from the 2-part model in that a latent variable is used to model whether the zeros arise from the occupancy or abundance components of the model.Note the difference in interpretation between the mixture and 2-part models: the mixture model estimates unconditional abundance; the 2-part model estimates conditional abundance, given presence.We fitted both model types and compared results for consistency.
A zero-inflated mixture model is formulated as follows.For response, y i , (Lambert 1992): where y represents densities, f (µ i ) is a distribution such as the binomial or Poisson and µ i is the mean.Co variates can be included in either or both components of the model: occupancy via the logit of the probability of a true zero, logit(θ i ), and abundance via a function of the mean, h(µ i ).The above formula is modified for the 2-part model such that µ i is the conditional mean and f (µ i ) is a zero-truncated distribution.The 2-part model used the proportion of surveys with a positive density in each cell as the Ά response for occupancy.Positive density values were the res ponse for abundance.The zero-inflated mixture model used whale presence or absence in a grid cell to model occupancy, and the density to model abundance.
Densities were log-transformed prior to analyses, with a constant of 0.1 added to the zero densities.The log-transformed densities were shifted to the right so that the smallest value (log (0.1)) was located at zero and the classification of a cell as having zero or nonzero density was maintained.
Spatial correlation was accounted for by including a CAR random effect (Besag et al. 1991) in one or both of the model components.The random effects b = (b 1 , b 2 , …, b N ), where each b j , j ⑀ 1, … N, is associated with a particular location, often account for heterogeneity from unmeasured spatially structured covariates.Under the CAR model, the random effects have univariate normal conditional distributions: where N is the number of grid cells, υ j = Σ i ⑀ δj w ij b i and δ j is a set of neighbours of grid cell j.We defined neighbours as adjacent grid cells and specified the weights w ij to be 1 for a neighbour and 0 otherwise.The Bayesian hierarchical models were run in Winbugs 1.4 (Spiegelhalter et al. 2003) using vague prior distributions.Fitting Bayesian hierarchical spatial models is relatively computationally intensive, especially for mixture models.To focus the Bayesian ana lyses, we conducted preliminary model exploration using generalized linear models and linear mixed effects models in a frequentist framework, with Akaike's Information Criteria (AICs) used to select a set of covariates for each of the occupancy and abundance components.These covariates were subsequently used in the Bayesian 2-part and mixture models.Details of the covariate selection modelling work are presented in Supplement 4 (www.int-res.com/articles/suppl/n029p211_supp.pdf).The co variate effects were given N(0, 10 000) priors, the intercept terms were given flat, uniform priors; the variance parameters associated with random effects were assigned inverse Gamma (0.5, 0.0005) priors or non-informative uniform priors.All simulations were run for an initial burn-in of 30 000 iterations.Additional samples were thinned to every 5 th sample until 20 000 samples had been obtained for estimating the posterior distribution.Point estimates for each covariate effect were based on the posterior mean, and 95% credibility intervals were based on the 2.5 th and 97.5 th percentile of the posterior distribution.

RESULTS
The seismic survey was conducted for 16 d from 17 June to 2 July 2010.There was considerable, albeit irregular, gray whale survey effort during this time.Foggy conditions frequently prevented surveys during the first half of the seismic activity (17−24 June), resulting in the majority of survey effort (~64%) occurring during the last half of the seismic activity.At least one complete shore-based survey of the 5 distribution stations was conducted on 10 of the 16 days.A total of 180 scans were conducted by distribution teams; 257 sightings of 315 gray whales were ob tained.Behaviour teams completed 110 scans that recorded 98 sightings of 126 gray whales.Most sightings (~80%) were of individual gray whales.Average daily counts of whale sightings per scan at each station showed an increasing trend in numbers throughout the monitoring (Muir et al. 2015a).Details of survey effort, number of gray whale sightings and individuals, and daily sighting maps are provided in Muir et al. (2015a).
The density surface consisted of 299 grid cells (Fig. 1).The percentage of non-zero density cells in the daily surfaces averaged 5.4 (SD: 3.03, range: 1.0 to 10.4).

Sound covariates
Only the cumulative SEL over the 3 h (cSEL3h) time window had zero values (~8%) because the 3 d and 7 d time windows spanned at least 1 sailing of a seismic line.Non-zero 3 h cumulative SEL (cSEL3h) had a mean of 152.7 dB re 1 µPa 2 -s (range: 54.8− 183.4).Three day and 7 d mean cSEL values across all time windows (cSEL3d, cSEL7d) were similar, with means of approximately 168 dB re 1 µPa 2 -s (range: 113.4−190.9).
The high ensonification level area (91 km 2 ) extended ~24 km along the portion of the study area closest to the seismic survey (Fig. 2).Moderate ensonification (87 km 2 ) occurred mainly to the north and south of the high ensonification area, and included a strip of cells in ~12 m depths between the high and low ensonification areas.Low ensonification (121 km 2 ) occurred mainly in the most southerly part of the study area and in very shallow depths; these areas are historically associated with low numbers of gray whales (Vladimirov et al. 2013).
Changes in cSEL over time were summarized well using FPCs for the 3 h time window.The first 3 FPC scores in the 3 h FPCA accounted for 83% of the survey-to-survey variation in sound pattern across grid cells.The first FPC (3hPC1) was positive and large when sound was low at the beginning, and high at the end of the 3 h (46% of total variation).The second FPC (3hPC2) was large when sound was high during the middle of the window (29% of total variation).The third component (3hPC3) was large when sound was low in the middle, and high at the beginning and end of the 3 h (8% of variation).
FPCs using an 8 h smooth for 3 d time windows did not seem to capture meaningful summaries of the data and were not used in any models.The first 3 FPCs using the 24 h smooth accounted for all of the survey-to-survey variation because there were three 24 h periods in the 3 d.While the FPCs were simply a transformation of the 24 h means, they also captured differences in the shape of the 24 h means in sound across the 3 d time windows.The first FPC (3d24h PC1) accounted for most of the variability in the three 24 h periods; it represented high sound at the beginning and end of the 3 d window.The second FPC (3d24hPC2) represented high sound at the beginning and low sound at the end.The third FPC score (3d24h PC3) represented a pattern of high sound in Day 2 (i.e. the previous 24−48 h), and low sound in Days 1 and 3 of the time window.
The large number of jumps in the longer 7 d windows did not produce easily interpretable FPCs and were not included in any models.

Occupancy
Both the 2-part and mixture models included a random effect for location (grid cell) and allowed the time effect (day category) to differ across grid cells.It was not feasible to include a CAR random effect in the occupancy models because ~95% of observed cell va lues were zero, resulting in perfect correlation among most adjacent cells.Estimates for the regression co efficients of the logistic occupancy models are presented as odds (Table 3).Increasing positive odds indicate a higher probability of occupancy.Posterior regres sion estimates were generally higher for the mixture model.The 2-part and mixture occupancy models both predicted the highest probability of occupancy in the moderate ensonification zone, followed by the high ensonification zone, and lastly the reference category low ensonification zone.However, there was considerable overlap in the high and moderate ensonification posterior distributions for both models (Fig. 3).
Both model types indicated similar temporal and spatial effects.The prob ability of cell occupancy in creased over time, reflecting the migration of gray whales into the study area.Probability of occupancy increased slightly with increasing northing (i.e. with more northerly grid cells), which is consistent with historically obser ved gray whale pref- Table 3. Two-part and mixture occupancy model posterior mean regression estimates as odds, with 95% credibility intervals shown in parentheses (2.5 th percentile, 97 th percentile).Odds of 1 indicate no effect.Models include fixed effects covariates and a random effect for grid cells.For ensonification categories, low ensonification represents the reference category that has odds of 1 Fig. 3. Two-part and mixture occupancy model posterior distributions for moderate and high ensonification regression estimates.The x-axis represents estimated log (odds) for these covariates in Table 3 erence for the northern part of the study area that is adjacent to the Piltun Lagoon.In addition, there was an interaction between easting (i.e.grid cell location in the east-west direction) and depth that reflected the complex bathymetry in the study area.The Bayesian model pre dictions are shown over time for the 2-part (Fig. 4) and mixture (Fig. 5) models.2).See Fig. 1 for explanation of stations

Abundance
Both models included a CAR spatial random effect and a random effect for day.A space-time effect was not included in the abundance models because there were few positive density cells during the first week of the seismic survey from which to estimate a spacetime interaction.The posterior mean regression coefficients for the fixed effects of the 2-part and mixture models are presented in Table 4. Higher sound levels in the previous 24 to 48 h combined with lower sound levels during the first and third day of the 3 d time window (3d24hPC3) were associated with a slight decrease in whale densities in both models.The 3d24hPC3 posterior distributions were nearly identical for the 2-part and mixture models (Fig. 6).Both models predicted slightly higher densities with increased northing (Fig. 6).As with the occupancy Fig. 5.As in Fig. 4, but for mixture model models, there was an interaction between easting and depth.Posterior means for the intercept terms and coefficients for the main effect of easting and the interaction between easting and depth were slightly smaller for the mixture model compared to the 2-part model.Predicted daily densities for the 2-part model are shown in Fig. 7.

Two-part and mixture model performance
Residual analysis had similar re sults for both 2-part and mixture models; only results for the 2-part model are presented.The occupancy model performed relatively well, correctly predicting whether or not a cell was occupied in the majority of cases, and there was a clear difference in the predicted probabilities of occu pan cy between occupied and unoccupied cells.Residuals for the abundan ce model did not appear to de viate from normality and no large outliers were observed.However, there was low correlation be tween observed and predicted values, indicating poor model prediction.

DISCUSSION
The primary objective of this study was to examine whether sound exposure from a 16 d seismic survey displaced gray whales from or within their feeding habitat.Gray whale distribution within our study area exhibits fine-scale variability that reflects inter alia the patchy nature of their benthic and epibenthic prey (Fadeev 2013, Vladimirov et al. 2013).Use of 1 km 2 density estimates for our analyses captured this spatial heterogeneity, and allowed estimation of more accurate sound covariates at grid cell centroids for association with 2 response variables, gray whale occupancy (presence/absence in density surface grid cells) and abundance (gray whale density within cells).Slightly decreased densities were as sociated with sound exposure in the previous 3 d that had a pattern of high sound in Day 2 and low sound in Days 1 and 3. Models also suggested highest occupancy in areas with moderate cumulative sound exposure levels.
The occupancy models showed higher odds of an occupied cell under moderate and high ensonification compared to low.In this context, it should be noted that even in the absence of additional noise, the geographic area covered by the low ensonification zone has very shallow waters and/or is in the far 222 Fig. 7. Two-part model predicted gray whale (Eschrichtius robustus) density on days with shore-based gray whale survey effort.Predictions were only made for occupied grid cells on a given day.Day 2 corresponds to 18 June 2010, which was the second day of seismic survey activity (fog on 17 June 2010 prevented shore-based observations of gray whales on that day).Model includes easting × depth, northing, 3d24hPC3 score (see Table 1 for definition), a random effect for grid cell and a conditional autoregressive (CAR) spatial random effect.See Fig. 1 for explanation of stations  south of the study area where there are typically low numbers of gray whales (Vladimirov et al. 2013).We believe this is a more likely explanation for a seemingly surprising result than an alternative explanation that changes in the whales' behaviour due to seismic activity caused them to spend more time on the water surface and hence be more visible.Gailey et al. (2016, this Theme Section) found no effects of sound exposure on 10 gray whale movement and 7 respiration response variables during the same seismic survey, although they cautioned that their analyses would not detect moderate to subtle effects.This issue warrants further consideration in future studies.
Our analysis required a complex modelling framework to address the trend of increasing numbers of gray whales as they migrated into the study area, considerable natural variability in gray whale distribution and abundance, and the large number of zero densities that remained even when multiple scan surveys were amalgamated into daily density surfaces.Use of zero-inflated models in ecological applications is widespread, as they are often needed for analyses of the abundance of rare species (Cunningham & Lindenmayer 2005).A second advantage of zeroinflated models is that they are decomposed into 2 parts: occupancy and abundance.This allows different covariates to be fitted for each component, which is often biologically relevant, as different covariates may relate to occupancy and abundance.
The development of zero-inflated models for continuous data was first proposed by Aitchison (1955) and Aitchison & Brown (1957), and a 2-part deltalognormal or delta-Gamma model (Fletcher et al. 2005) has been used in ecological applications (e.g.Pennington 1983, Stefánsson 1996).We implemented a 2-part delta-lognormal model and the analogous zero-inflated mixture model in a Bayesian spatiotemporal framework to assess effects of characteristics of sound exposure on gray whale occupancy and abundance.The key difference between these models is that the mixture model separates 'true' zeros (i.e.arising from unsuitable habitat) from 'false' zeros (due to temporary absence of whales or a failure to detect whales that were present during a scan) into the occupancy and abundance components, respectively, of the model (Martin et. al. 2005).In contrast, the 2-part model does not distinguish different types of zeros and includes all zeros in the occupancy component, with an abundance component consisting only of positive densities.We therefore fitted both 2part and mixture models so that results from both model types could be evaluated for consistency and if found to be so, provide more confidence in infer-ences.Results were similar, although the mixture model tended to predict slightly higher probabilities in the occupancy component and lower probabilities in the abundance component that reflected the separation of true and false zeros in these components.However, the 2-part models were simpler to fit and converged more easily.Both model types revealed similar effects of sound and a strong association of both gray whale occupancy and abundance with geographic location and depth.The probability of occupancy increased over time, reflecting the trend of increasing numbers of gray whales as they mig rated into the Piltun feeding area.
Our findings suggest a possible avoidance re sponse by gray whales to higher sound exposure levels after a prolonged period of disturbance.However, we recognize that effects on gray whale densities need to be interpreted with caution, given the limited positive abundance data available with which to make inferences.Although we conducted a large number of scans at the observation stations, there were few gray whales in the study area during the monitoring, and hence most sampled grid cells had zero density estimates.Consequently, there was limited information in the collected data for estimating the relationships between occupancy and/or density and sound covariates.The abundance models, in particular, had low correlation between observed and predicted values, indicating poor model prediction.It is also possible that whales were responding to other factors, including changing prey availability, al though the time interval was short (16 d).Data on prey distribution and abundance were not available and thus this potentially important covariate could not be included in our models.Intra-seasonal shifts in gray whale distribution hypothesized to be associated with varying prey supply are commonly seen within our study area (Fadeev 2013, Vladimirov et al. 2013).Similarly, Dunham & Duffus (2001) found that foraging gray whales off British Columbia, Canada, moved throughout their feeding area in response to changes in prey supply.As had originally been plan ned as part of the monitoring programme (Bröker et al. 2015), pre-and post-seismic data might have allowed effects from sound exposure and prey availability to be separated.Although monitoring was conducted before and after the seismic survey with the intent of including a before/after impact component in the analysis, limited data were obtained.Very low numbers of gray whales were present during pre-seismic monitoring that coincided with the be ginning of gray whale migration into the feeding area, resulting in density surfaces with extremely high numbers of zeroes.Foggy conditions severely limited survey effort during most of the planned 2 wk post-seismic monitoring when overall whale numbers were much higher.We therefore could not de termine conclusively whether observed changes in occupancy and abundance were due to sound exposure levels or to natural variation in gray whale distribution.
Other studies have also related sound exposure over a 3 d period to observed shifts in baleen whale distribution.Yazvenko et al. (2007) attributed an observed redistribution of an estimated 5 to 10 whales in their study's feeding area to a higher 3 d cumulative sound exposure from a seismic survey.Castellote et al. (2012) similarly found that fin whales Balaenoptera physalus, as indicated by number and received levels of song notes, moved away from an airgun array and out of the study's detection range over a 3 d period.This displacement lasted until approximately 14 d after seismic activity had ceased.Our results suggest that sound exposure over the preceding 3 d can affect gray whale distribution and abundance, and a lagged effect of higher sound levels on the previous day may be important.Further work to analyse effects of cumulative sound levels within this specific time period would be informative.
Our study was distinctive in that we modelled time series of 5 min cumulative sound exposure bins at each density surface grid cell to estimate sound covariates that allowed testing of the effects of both total magnitude and characteristics (encroachment pattern, proportion, median levels and variability of 'on' sound) of sound exposure from the airguns over time windows of 3 h, 3 d and 7 d preceding each gray whale survey.Additionally, we used FPCA to summarize the patterns in time series of cumulative sound exposure levels at cell locations into a small number of variables that were used as covariates in our models.FPCA worked well over 3 h, but had limited success for 3 d and 7 d time windows that had several on/off jumps in sound exposure because the seismic vessel sailed multiple lines during these time periods.To address this issue, we used FPCA on 24 h averages within 3 d windows to model broad changes in the pattern of sound exposure levels from day to day.FPCA is a relatively new technique, and methods to handle jumps in time series are the topic of ongoing research.FPCA has been successfully used to relate temporal variation in river flow to ecological responses such as fish abundance (Ainsworth et al. 2011, Stewart-Koster et al. 2014), and thus could be useful for assessing effects of continuous sound exposure (e.g.offshore construction, vessel traffic) on marine mammal occurrence or abundance.
We recommend further exploration of the analytical approach developed in this paper for future seismic surveys in the area.Under certain assumptions it may be possible to combine data collected in different studies and conduct a meta-analysis.The in creased numbers of non-zero densities in the pooled dataset will provide more information to estimate the relationships between sound covariates and re sponses of occupancy and abundance.However, for future studies, we caution that mitigation to reduce sound exposure remains a priority over increasing non-zero samples by conducting a seismic survey when the Piltun area has a higher abundance of gray whales.

Fig. 1 .
Fig. 1.Study area, 4-D seismic survey area (green polygon with black seismic lines) and grid cell coverage for Distribution Stns 9-13 (black triangles) and behaviour Seismic North and Seismic South (purple circles) shore-based stations that conducted scans of gray whales Eschrichtius robustus that were used to estimate densities.Estimated 20 m and 50 m bathymetry contours are shown as dotted blue and solid blue lines respectively

Fig. 4 .
Fig. 4. Two-part model predicted probability of grid cell occupancy by a gray whale Eschrichtius robustus for each day category.Model includes easting × depth, northing, day category, ensonification, random effects for grid cell and day category coefficient.Day 1 is 17 June 2010, the first day of seismic activity.Day number groups shown above the panels correspond to the 3 day categories used in the model (Table2).See Fig.1for explanation of stations

Table 4 .
Two-part and mixture abundance model posterior mean regression estimates for fixed effects, with 95% credibility intervals shown in parentheses (2.5 th percentile, 97 th percentile).Models include fixed effects covariates, random effect for grid cells and a conditional autoregressive (CAR) model of spatial random effect.See Table1for definition of 3d24hPC3