The population density of Cuvier's beaked whales is estimated acoustically with drifting near-surface hydrophone recorders in the Catalina Basin. Three empirical approaches (trial-based, distance-sampling, and spatially explicit capture-recapture) are used to estimate the probability of detecting the echolocation pulses as a function of range. These detection functions are used with two point-transect methods (snapshot and dive-cue) to estimate density. Measurement errors result in a small range of density estimates (3.9–5.4 whales per 1000 km2). Use of multiple approaches and methods allows comparison of the required information and assumptions of each. The distance-sampling approach with snapshot-based density estimates has the most stringent assumptions but would be the easiest to implement for large scale surveys of beaked whale density. Alternative approaches to estimating detection functions help validate this approach. The dive cue method of density estimation has promise, but additional work is needed to understand the potential bias caused by animal movement during a dive. Empirical methods are a viable alternative to the theoretical acoustic modeling approaches that have been used previously to estimate beaked whale density.

There is a growing interest in using passive acoustic methods to estimate the population density and abundance of cetaceans (Marques et al., 2013; Heinemann et al., 2016). Previously, most cetacean abundance estimates were based on visual sighting surveys (Dawson et al., 2008) or capture-recapture methods with photographic identification of individuals (Urian et al., 2015). However, some species of cetaceans, such as beaked whales, are difficult to see or photograph because they spend very little time at the surface and are relatively inconspicuous (Barlow and Gisiner, 2006). Because beaked whales make relatively loud echolocation pulses (224 dBpp re 1 μPa; Gassmann et al., 2015) when foraging and because these pulses appear to have species-specific characteristics (Baumann-Pickering et al., 2013), passive acoustic abundance estimation may be especially effective for these species. In some Navy test facilities with large arrays of permanent, bottom-mounted hydrophones, it is possible to estimate abundance by completely enumerating the number of beaked whale foraging dives and adjusting for the dive rate and group size (Moretti et al., 2010; Marques et al., 2019). However, for most areas, fewer hydrophones are available, and a sampling approach must be taken. Various approaches have been proposed (Marques et al., 2013), most of which are based on point-transect sampling.

A key component of point-transect sampling for animal density or abundance is estimating the probability of detecting an animal or group of animals as a function of their distance from a point (e.g., a single hydrophone). This function, referred to here as a detection function, can be estimated by a variety of approaches. One approach is based on theoretical models that use estimated sound levels at the source, propagation losses with range, and the required signal-to-noise level at the hydrophone for an acoustic signal to be recognized. This approach is more complicated for higher-frequency sounds that are highly directional and have a broad bandwidth (Au, 1993; Zimmer et al., 2008; Ainslie, 2013). For such highly directional sounds, such as beaked whale echolocation pulses, a detailed simulation model of animal behavior is also needed to predict the detection range (Zimmer et al., 2008; Küsel et al., 2011; Hildebrand et al., 2015). However, these acoustic models are typically based on a large number of untested assumptions and unquantified sources of uncertainty regarding individual and group foraging behaviors. Empirical methods for estimating detection functions include distance-sampling, mark-recapture distance-sampling, trial-based, and spatially explicit capture-recapture (SECR) approaches (Marques et al., 2013). Empirically estimated detection functions are likely to be more robust (Marques et al., 2013) but require a survey design that collects information on detection distances.

Conventional distance-sampling uses only the empirical distribution of ranges to actual detections and does not require knowledge of the range to missed detections (Buckland et al., 2015). A detection function is estimated from the observed decrease in the number of detections with range (relative to the expected number of detections if range did not affect detectability). However, the empirical estimation of range to an acoustic detection is not a trivial problem, especially when using a single hydrophone. In some special cases, range from a single hydrophone to an acoustic source can be estimated using surface and bottom reflections (McDonald and Fox, 1999) or waveform dispersion caused by waveguide propagation (Wiggins et al., 2004; Marques et al., 2011). Additional methods are available for two hydrophones in a vertical array (Mathias et al., 2013; Barlow and Griffiths, 2017). Distance-sampling requires an independent measurement of the probability of detection at zero distance (or the assumption that detection is certain) and typically assumes that detection probability decreases monotonically with distance.

The trial-based approach uses both detections and missed detections to estimate the detection function. In a passive acoustic application, this approach requires that the location of the animal be estimated both when signals are being received and when they are not being received on a given instrument. The information on the range to missed detections can come from other instruments which did detect the animal (such as an acoustic recording tag on the animal; Marques et al., 2009) or from a modeled track of animal movement that estimates locations over a continuous time period from intermittent acoustic detections. The trial-based method used with nonparametric smoothing functions does not assume that detection is certain at zero distance or declines monotonically.

The SECR approach estimates a detection function using multiple receivers distributed in a known spatial arrangement. When the same call is detected on multiple receivers, each detection is essentially a “recapture” event, and the source location of the detected call can be estimated from the pattern of recaptures (Marques et al., 2012; Stevenson et al., 2015). SECR methods do not (necessarily) assume that detection is certain at zero distance; they do typically assume a monotonic decreasing detection function because of the parametric detection functions used, although non-monotonic functions may be possible in theory.

The effective detection radius (EDR) is a useful summary statistic for the detection function. The EDR is the distance at which the number of missed detections at closer distances equals the number of detections made at greater distance (within a defined maximum truncation distance). This can be readily converted into the effective area surveyed (π×EDR2), which is a key component of density estimation (see Sec. II, Methods). For the SECR approach, this area depends on the spatial distribution of instruments.

Much recent research has addressed the range at which echolocation pulses of beaked whales can be detected using acoustic modeling and empirical approaches. Zimmer et al. (2008) used a sound propagation theory and a simulation model of animal behavior to estimate the probability of detecting a Cuvier's beaked whale (Ziphius cavirostris) as a function of range within a finite time period (taken to be the acoustically active phase of one foraging dive). Sound propagation theory and a simulation model of animal behavior were used to estimate the probability of detecting individual echolocation pulses for Blainville's beaked whales (Küsel et al., 2011), Cuvier's beaked whales, and Gervais' beaked whales (Mesoplodon europaeus; Hildebrand et al., 2015). Hildebrand et al. (2015) also used an acoustic model to estimate the probability of detecting Cuvier's and Gervais' beaked whale groups within a five-minute time window. Marques et al. (2009) used a trial-based approach to empirically estimate the probability of detecting individual echolocation pulses for Blainville's beaked whales using pulses detected on an acoustic recording tag to set up trials for detection or non-detection on bottom-mounted hydrophones at varying distances. A fundamental difference in the unit of measure (detection of a dive vs a finite time window vs an individual echolocation pulse) prevents a detailed comparison of the EDR from these papers, but Zimmer et al. (2008), Küsel et al. (2011), and Hildebrand et al. (2015) all indicate a maximum likely detection range of slightly greater than 4 km, whereas Marques et al. (2009) extend that to a little over 6 km.

In this paper, three approaches (distance-sampling, trial-based, and SECR) are used to empirically estimate the acoustic detection function for Cuvier's beaked whales on drifting near-surface hydrophones. Acoustic detection data are taken from a study that used a nested array of drifting hydrophone recorders to track the three-dimensional (3D) diving behavior of this species (Barlow et al., 2018). Estimated detection functions are used to estimate the effective survey area of the drifting instruments, and that area is used to estimate the density of Cuvier's beaked whales in a small study area.

We use acoustic detections of echolocation pulses from Cuvier's beaked whales that were recorded as part of a beaked whale acoustic tracking study (Barlow et al., 2018). That study was conducted in the Catalina Basin off southern California in July and August 2016 and used drifting arrays of hydrophones to localize and track foraging whales based on their echolocation pulses. Detailed field methods are given by Barlow et al. (2018); here, we present a brief summary.

Acoustic detections are from four Drifting Acoustic Spar Buoy Recorders (DASBRs). Each consists of a vertical array of two hydrophones at ∼105-m and 115-m depths (separated by 10 m) and a Wildlife Acoustics (Maynard, MA, USA) SM3M digital recorder with a 256 kHz sampling rate. The instruments drifted together with prevailing currents and were repositioned daily with an initial separation of approximately 1 km. The location of each instrument was recorded approximately every 30 min with a Global Positioning System (GPS)-based geo-location device. Drift paths were generally in a northwesterly direction during this study (Fig. 1). Eleven drifts (a single deployment of one instrument) were completed with each of the four recorders.

FIG. 1.

(Color online) Locations of DASBR drifts in the Catalina Basin (black lines). Four DASBRs were typically deployed in a rectangular configuration separated by ∼900 m (as exemplified by black circles showing deployment locations for one drift) and drifted northwest. Locations of Cuvier's beaked whales are illustrated as white squares with black outlines. Bathymetry contours are every 200 m from 200 to 1200 m (light gray lines) with the 1000 m contour emphasized with the darker, thicker gray line. Bathymetry data are from National Oceanic and Atmospheric Administration's (NOAA's) National Centers for Environmental Information (Amante and Eakins, 2009). The inset shows the location of the detailed map within the Southern California Bight.

FIG. 1.

(Color online) Locations of DASBR drifts in the Catalina Basin (black lines). Four DASBRs were typically deployed in a rectangular configuration separated by ∼900 m (as exemplified by black circles showing deployment locations for one drift) and drifted northwest. Locations of Cuvier's beaked whales are illustrated as white squares with black outlines. Bathymetry contours are every 200 m from 200 to 1200 m (light gray lines) with the 1000 m contour emphasized with the darker, thicker gray line. Bathymetry data are from National Oceanic and Atmospheric Administration's (NOAA's) National Centers for Environmental Information (Amante and Eakins, 2009). The inset shows the location of the detailed map within the Southern California Bight.

Close modal

Cuvier's beaked whales were identified from their characteristic echolocation pulses (Baumann-Pickering et al., 2013). Digitally recorded .wav files were analyzed using PAMGuard (Beta v1.15.03) software (Gillespie et al., 2009), which identified impulsive sounds (clicks and pulses) using an energy detector and classified those sounds based on peak frequencies and, for some, the presence of a frequency sweep within an echolocation pulse (Keating and Barlow, 2013). PAMGuard also estimates a vertical bearing angle to the sound source based on the time-difference-of-arrival for the signal on the two hydrophones in the vertical array using a cross-correlation algorithm. Potential echolocation pulses from beaked whales were reviewed by an experienced analyst (E.T.G.), who used PAMGuard Viewer software to determine whether signals are from Cuvier's beaked whales, based on their frequency characteristics and contextual cues. Echolocation pulses of Cuvier's beaked whales characteristically have a peak frequency of 32–40 kHz and secondary energy peaks at 18–19 kHz and 22–24 kHz. As with other beaked whales, the echolocation pulses of Cuvier's beaked whales have a frequency upsweep. Important contextual cues include the direction to the sound source (always below the hydrophone pair), the consistency of that direction (changing only slowly over a two-minute sound file), and the inter-pulse interval (typically 0.33–0.50 s). Because not all echolocation pulses can be classified unequivocally, three clicks within a one-minute sound file are required to qualify as a positive detection of a beaked whale. If a series of potential beaked whale echolocation pulses differs from those of Cuvier's beaked whales in any of the diagnostic criteria, it was excluded from subsequent analyses. Our multiple-criteria approach is assumed to eliminate any false positive detections; however, false negatives (missed detections) are likely and are the motivation for estimating detection functions.

Two methods are used for acoustic density estimation, both based on a point-transect framework (Marques et al., 2013; Buckland et al., 2015), using groups as the sampling unit. The first method uses a one-minute “snapshot” sample during which a group of Cuvier's beaked whales can be either detected or not. In point-transect surveys, the choice of snapshot duration should be long enough to have a high probability of detecting close groups (as we will show to be true for one-minute snapshots) but should be as short as possible to reduce the bias caused by animal movement and the likelihood of detecting multiple groups (Buckland et al., 2015). Because fewer than 2% of the snapshots include whale detections, we are confident in assuming that a one-minute snapshot detection includes only one group of beaked whales. The density of individuals, D, is estimated as

(1)

where n is the number of snapshot time periods with detections of Cuvier's beaked whales, k is the total number of snapshot time periods, ν̂ is the estimated effective area surveyed (which includes the effect of the detection function), λ̂ is the estimated probability that a group will be available for detection (echo-locating) within a one-minute snapshot, and ŝ is the estimated population mean group size. The mean group size (1.9, coefficient of variation, CV = 0.07) is estimated as the mean of the observed group sizes from 63 Cuvier's beaked whale sightings on 7 prior visual sighting surveys in the California Current (Barlow, 2016). The instantaneous proportion of time that an individual Cuvier's beaked whale echo-locates (0.199, CV = 0.04) is taken from a tagging study in similar habitat in southern California (Barlow et al., 2020). The probability, λ̂, of a Cuvier's beaked whale group echo-locating during a one-minute time window is approximately one percentage point higher than this proportion of time echo-locating (Barlow et al., 2013). Individuals in a group are assumed to begin and end echolocation synchronously, thus, from this instantaneous proportion, availability for a one-minute snapshot, λ̂, is estimated as 0.209 (CV = 0.04). Three methods are used to estimate the probability of detection as a function of horizontal range (Sec. II D) and the effective detection area ν̂ (Sec. II F).

The second density estimation method is cue counting using foraging dives as cues. The density of individuals, D, is estimated as

(2)

where nd is the number of detected dives of Cuvier's beaked whales during a time period T hours, ŝ is the population mean group size, ν̂ is the estimated effective area surveyed, and R̂d is the estimated dive rate (dives hr−1). This method is based on that of Moretti et al. (2010) except that the Moretti study used data from a Navy testing range and was able to assume that all dives within a known area were detected, whereas here detection of group dives is not certain and so the effective area surveyed is estimated. The group size is again taken as 1.9 (CV = 0.07) from visual sighting surveys. Again, at least three beaked-whale-like echolocation signals detected within a one-minute period are required to qualify as a detection. In this region, individual Cuvier's beaked whale foraging dives have a mean duration of 67.4 min (standard deviation, s.d. = 6.9) and a mean period between foraging dives of 102.3 min (s.d. = 30.8) (Schorr et al., 2014). This species typically does not produce echolocation clicks during the descent or ascent phases of these dives, therefore, the period of echolocation is considerably shorter than the length of a foraging dive. A dive is counted in nd if the time since the detection of a previous dive was greater than 50 min. The dive rate (0.319 dives hr−1, CV = 0.13) is taken from a southern California tagging study (Barlow et al., 2020). Each detected dive is assumed to contain only one group, and individuals within a group are assumed to be diving synchronously. Two dive detection functions are estimated (Sec. II E) to bracket the likely range of EDR values, and both are used to estimate ν̂ (Sec. II F). Animal movement during a dive is likely to bias density estimates by this method (see Sec. IV, Discussion).

The precision of both density estimation methods is estimated with the delta method. The square of the CV of the density is estimated as the sum of the squared CVs of each multiplicative component in Eqs. (1) and (2).

Three approaches are used to estimate the probability of detection as a function of range (trial-based, modified distance-sampling, and SECR). The trial-based approach only used data from 10 dives of solitary individuals that could be tracked by Barlow et al. (2018), whereas the other 2 approaches used data from all 29 dives detected in that study, some of which included multiple individuals. For all approaches, the EDR is estimated from detection functions using the identity

(3)

where r is the horizontal range in meters, g(r) is the probability of detection at range r, and w is a range beyond which detection probability is essentially zero (Buckland et al., 2015). Given the knowledge of the beaked whale acoustic detection (Zimmer et al., 2008), the chosen value of w was 6000 m. In practice, the integration was approximated by a sum taken over integer ranges. The estimated effective surveyed area, ν̂ Eqs. (1) and (2), is calculated as the area of a circle with a radius of the EDR:

(4)

1. Trial-based detection probability

The trial-based approach uses data on whether an animal at a given horizontal range and orientation was detected or not within a one-minute snapshot on each of the four DASBRs (Fig. 2). Barlow et al. (2018) developed a Bayesian state-space model to estimate the locations of individual Cuvier's beaked whales at 1-minute intervals during 10 dive segments of 9–32-minute duration using a subset of the data used in this study (see the supplementary material in Barlow et al., 2018). These estimated locations were used to compute the horizontal range and relative heading from each DASBR. The relative heading is informative because beaked whale echolocation pulses are highly directional and, hence, the detection range depends on an animal's orientation relative to the hydrophone. Although an animal's precise orientation is not known, it is likely to be correlated with its direction of travel and, therefore, the horizontal orientation is estimated from the time series of tracking locations and converted to a heading in degrees relative to each DASBR (with 0° being directly toward the DASBR and 180° directly away from the DASBR).

FIG. 2.

(Color online) Acoustic detections (red circles) and non-detections (black circles) on a DASBR during a one-minute time window as a function of horizontal distance from that DASBR and heading of travel relative to that DASBR. Data are from an acoustic tracking study (Barlow et al., 2018).

FIG. 2.

(Color online) Acoustic detections (red circles) and non-detections (black circles) on a DASBR during a one-minute time window as a function of horizontal distance from that DASBR and heading of travel relative to that DASBR. Data are from an acoustic tracking study (Barlow et al., 2018).

Close modal

A generalized additive model (GAM) is fit to estimate the per-minute detection probability as a function of the range and relative heading. Using GAMs avoids the need to specify a parametric form for the function; instead, it is modeled using spline functions, which are smooth but flexible (Wood, 2017). The model does not assume that the probability of detection at zero range is one, nor does it assume that the detection probability decreases monotonically with range. The response variable is detection/non-detection per minute, modeled as a Bernoulli trial with a logit link function; the two explanatory variables (horizontal range and relative heading) are included as univariate thin plate regression splines. Fitting is via maximum likelihood with the amount of smoothness of the splines being selected as part of the maximization. In doing so, it is important to recognize that the binary trials on each DASBR in each minute are not statistically independent because each minute of tracking generates four trials (one on each DASBR) and trials within tracks are based on the same diving group. Hence, independent random effects were included in the model for a minute within track and for DASBR. Fitting was performed using the R (R Core Team, 2017) function gamm4 (version 0.2–5; Wood and Scheipl, 2017; see supplementary material 1 for R code1).

To derive the detection probability as a function of the horizontal distance alone, the relative heading was “integrated out” by assuming a uniform distribution of headings between 0° and 180° and averaging the detection probabilities for all angles at a given range. The EDR was then calculated using Eq. (3). For increased robustness, the variance was estimated using a jackknife procedure (see Sec. II F).

2. Distance-sampling detection probability

A variety of parametric and semi-parametric models can be used to model a detection function in distance-sampling (see, e.g., Buckland et al., 2015), and these are typically fitted to the observed distance data using the maximum likelihood. In our case, with acoustic data from drifting buoys, detection distances are not directly observed. Instead, our observations are of detection angles (in a vertical orientation). However, by combining these angles with separate information about the distribution of echolocation depths, it is possible to make inferences about the distribution of the detection distance and, hence, detectability. The mean depth of echolocation (976 m, s.d. = 112 m) is taken from the ten dive segments for which tracks were available (Barlow et al., 2018), and the distribution of depths was assumed to be log-normal. Rather than the maximum likelihood, the maximum simulated likelihood is used to find the detection function that provides the best fit to the observed angles (see supplementary material 2 for details1). The maximum simulated likelihood was first developed in the field of economics as an approach for fitting complex likelihoods (Arias and Cox, 1999). In the present application, the full likelihood involves integration over the depth distribution, and the simulation approach avoids this integration step [although full likelihood methods have been developed for the half-normal model (Thomas and Matias, 2014)—see Sec. IV, Discussion]. Fitting was done using custom-written code in R (see supplementary material 3 for R code1).

Half-normal (Buckland et al., 2015) and compound half-normal (Efford and Dawson, 2009) detection functions (Table I) are fit to the observed distribution of detection angles. Unlike conventional distance-sampling, both are modeled as functions of the slant (i.e., direct) range to detection rather than the horizontal range. A simulation is used to re-express these detection probabilities as functions of the horizontal range for estimating the EDR (see supplementary materials 2 and 31). Detection angles are truncated at 81.8° (corresponding to the horizontal truncation distance of 6000 m at the mean echolocation depth of 976 m). The best model is selected as the model with the lowest AIC (Akaike's information criteria) value.

TABLE I.

Maximum simulated likelihood fits of the detection model for the distance-sampling approach. Detection probability is a function of the slant range, r. In the formulas, σ and z are estimated parameters. The effective detection radii are given for the maximum simulated likelihood fits to each model. Relative model fits are indicated by Akaike's information criteria (AIC) and the probabilities are indicated from a Kolmogorov-Smirnov test [Pr(K-S)] of the differences between the observed and simulated distributions of detection angles.

Detection modelFormulaEDR (m)AICPr(K-S)
Half-normal er2/2σ2 2254 4099.7 0.0003 
Compound half-normal 1(1e(r2/2σ2))z 3000 4072.1 0.529 
Detection modelFormulaEDR (m)AICPr(K-S)
Half-normal er2/2σ2 2254 4099.7 0.0003 
Compound half-normal 1(1e(r2/2σ2))z 3000 4072.1 0.529 

3. SECR detection probability

Traditionally, SECR uses data from a set of physical traps that are deployed to catch animals over a set of discrete time periods called sampling occasions. Captured animals are individually identifiable or marked on first capture to make them so. The data are summarized in the form of capture histories, which records the history of which traps each detected animal was caught in and on which occasions. The number and pattern of detection locations can be used to estimate the detection function and effective area surveyed and, hence, density (see Borchers, 2012, for an introductory overview of SECR). When applying SECR in the context of passive acoustics, sounds act as animals and hydrophones act as traps—the number of individual sounds detected and the pattern of detection locations are used to estimate an acoustic detection function, effective area, and density (Efford et al., 2009b). Unlike with physical traps, a single sound can be detected on multiple hydrophones in the same time window, and so only one sampling occasion is needed.

SECR surveys can be divided into sessions (groups of sampling occasions separated by time or space), and the effective survey area and density are calculated for each session. This allows for comparison across sessions (Efford et al., 2009b; Royle et al., 2013). Conversely, sessions can be collapsed together to look at the mean density for all sessions (e.g., Marques et al., 2012). Sessions can also be useful because they allow for changes in trap or recorder configurations in time (Efford, 2019a) as occurred with the DASBR array.

Similar to the trial-based method, an upper limit to the feasible detection distances must be specified in the form of a buffer around the sensors—this defines the area, or mask of grid points, over which the likelihood is integrated (Borchers and Efford, 2008). Density estimates should not change with increased buffer size if the buffer is sufficiently large (Royle et al., 2013).

Detection probability as a function of the horizontal range is modeled in the R package secr (v3.2.1; Efford, 2019b; see supplementary material 4 for R code1). Each minute of recording with at least three beaked whale clicks on at least one instrument (n = 275) makes up a single capture history with the binary presence or absence of clicks on each of the four recorders. To account for movement of the array through time, each single capture history is treated as a separate session with a unique array configuration. The half-normal and compound half-normal detection functions (defined in Table I) are fit by maximum likelihood using the more stable conditional likelihood option in secr (Borchers and Efford, 2008).

A secr buffer size of 8 km is used as a generous buffer beyond estimated detection ranges found with the distance-sampling and trial-based methods. AIC is used to identify the best model. Both detection function parameters [including the probability of detection at horizontal distance zero, g(0)] are estimated by the secr program. SECR estimates the effective survey area for the entire array rather than a single instrument as is the case with distance-sampling and the trial-based approaches. The array configuration may not lead to an exactly circular effective survey area, therefore, an estimate of the EDR is not typically calculated. In the present case, for comparison, a “pseudo” effective detection range (pEDR) is calculated for a single instrument of the array using the mean of the jackknife estimates for the SECR parameters g(0), σ, and z (see Table I for definitions) to define g(r) in Eq. (3).

Our dive-based cue-count approach requires an estimate of the effective area surveyed or, equivalently, the EDR. The EDR is expected to be greater for a dive than for a one-minute snapshot because there are many more opportunities (minutes) to detect a dive. Also, animal orientation relative to the hydrophone is likely to change during a dive (Barlow et al., 2018), which would also increase the probability that a dive will be detected. Based on measurements taken from Cuvier's beaked whales, tagged off southern California (Barlow et al., 2020), the average duration of the foraging portion of a dive (when echolocation pulses are made) can be estimated from the mean foraging dive duration (65.5 min) minus the descent time to the depth where echolocation begins (462 m at 1.45 m s−1 equals 5.3 min) and minus the ascent time from the depth where echolocation ends (881 m at 1.45 m s−1 equals 21.0 min). Therefore, this echolocation period is approximately 39 min for an average dive. If the probability of detection during each minute were assumed to be independent, the probability of detecting a dive at a given range could be estimated as the converse of the probability that echolocation pulses would not be detected at that range for 39 consecutive minutes. If 1 - g(r) is the probability of not detecting an echolocating whale during a 1-minute period at range r, the probability of detecting a 39-minute dive at that range would be 1 - (1 - g(r))39. Because the assumption of independence is not likely to be met, the true dive detection probability is likely to be less than that predicted by this formula. In applying our cue-count approach, this 39-minute detection function and the 1-minute detection function are used to bracket the feasible range of dive detection functions. The 1-minute detection function from the distance-sampling approach with a compound half-normal model is used to estimate the 39-minute detection function.

Calculation of the effective area surveyed, ν̂, depends on the method used to estimate the detection function. For the trial-based and distance-sampling approaches, the effective survey area is measured as the area around each instrument [Eq. (4)].

For the SECR approach, the effective survey area is the total area for all four instruments and depends on the spatial distribution of the instruments. This area is estimated within the program secr.

To calculate the CV of all estimates of the effective survey area, a jackknife approach (Efron, 1982) was used with each of the 11 drifts as replicates (for the distance-sampling and SECR approaches) or the 10 tracked individuals as replicates (for the trial-based approach).

1. Trial-based approach

The tracking methods used by Barlow et al. (2018) resulted in 224 estimated locations in 1-minute intervals over the 10 tracks. Because the tracks were continuous and positions were interpolated based on a movement model even when no echolocation pulses were received, these locations included both one-minute periods with and without recorded pulses. The actual number of detection-positive minutes was 135, 146, 117, and 167 for the four DASBRs. These 565 total detection-positive minutes, thus, represent 63% of the 896 potential periods for the 4 DASBRs and 224 localizations used in this study. The horizontal distances from the localized positions to a DASBR ranged from 0.53 to 4.45 km; however, no detections were made at distances over 4.15 km. Detections and non-detections are shown in Fig. 2 as functions of the horizontal distance and direction of travel relative to a DASBR. These results show that most tracked groups had a net travel heading that was toward the DASBRs (<90°), likely because such groups were more trackable (a potential bias in the trial-based method—see Sec. IV, Discussion). However, even when the net travel heading was away from the DASBR (>90°), beaked whales within 1 km were very likely to be detected.

These results are reflected in the GAM fit, where estimated detection probability generally declines with an increasing range and increasing orientation angle (Fig. 3), although the decrease with the range is not monotonic. The estimated detection probability, averaged over all travel directions (Fig. 3), is 1.0 at zero horizontal range and decreases non-monotonically at greater ranges; the estimated detection probability is essentially 0.0 at ranges greater than 5 km. The estimate of the EDR is 2.79 km (jackknife CV = 0.19).

FIG. 3.

Trial-based estimates of the detection probability as a function of the horizontal range and relative heading from a GAM fit to detection and non-detection data (Fig. 2). Solid lines (gray scale) represent probabilities of detection in a one-minute snapshot for orientations of 0°, 30°, 60°, 90°, 120°, 150°, and 180° (respectively, from black to light gray). The dashed line represents the one-minute detection probability averaged over all headings.

FIG. 3.

Trial-based estimates of the detection probability as a function of the horizontal range and relative heading from a GAM fit to detection and non-detection data (Fig. 2). Solid lines (gray scale) represent probabilities of detection in a one-minute snapshot for orientations of 0°, 30°, 60°, 90°, 120°, 150°, and 180° (respectively, from black to light gray). The dashed line represents the one-minute detection probability averaged over all headings.

Close modal

One minor issue noted during the jackknife procedure was that the estimated average detection probability at 6 km was not close to zero [ĝ(w)=0.079] in one of the ten jackknife replicates. Given that our limit of integration for estimating the EDR is 6 km, this will cause a slight underestimation of the EDR for that replicate and, thus, a slight underestimation of the jackknife CV. [This was not an issue in the other nine jackknife replicates: mean ĝ(w) over these nine was <0.001.]

2. Distance-sampling approach

The distribution of mean detection angles for each one-minute snapshot (Fig. 4) shows a peak at 70°, which corresponds to a horizontal distance of 2.4 km for an animal at the mean depth of 967 m. The truncation angle (81.8°, corresponding to a horizontal distance of 6.0 km at this mean depth) excludes only 0.7% of the observed detection angles. The compound half-normal is better than the half-normal both in terms of having a lower AIC and higher Kolmogorov-Smirnov (K-S) goodness of fit p-value (Table I; see supplementary material 5 for cumulative distribution plots1). The compound half-normal shows a high probability of detection (nearly certain) out to approximately 2 km, whereas the simple half-normal shows an initial detection probability of approximately 0.88 at zero horizontal distance (Fig. 5). Both models show a low probability of detection at 4 km (<0.1) and a very low probability of detection at 5 km (<0.05). The compound half-normal estimate of the EDR is 3.00 km (jackknife CV = 0.15).

FIG. 4.

The observed distribution of beaked whale acoustic detection angles for one-minute snapshots. Zero corresponds to straight down. The dotted vertical line indicates the truncation angle of 81.7°; greater angles were not used to fit the detection functions for the conventional distance-sampling method.

FIG. 4.

The observed distribution of beaked whale acoustic detection angles for one-minute snapshots. Zero corresponds to straight down. The dotted vertical line indicates the truncation angle of 81.7°; greater angles were not used to fit the detection functions for the conventional distance-sampling method.

Close modal
FIG. 5.

Distance-sampling estimates of the detection probability as a function of the horizontal range for 1-minute snapshots and for a 39-minute dive. Detection models are half-normal (dashed) and compound half-normal (solid). The dotted line shows the expected detection function for a dive based on 39 independent 1-minute observations with the compound half-normal detection function.

FIG. 5.

Distance-sampling estimates of the detection probability as a function of the horizontal range for 1-minute snapshots and for a 39-minute dive. Detection models are half-normal (dashed) and compound half-normal (solid). The dotted line shows the expected detection function for a dive based on 39 independent 1-minute observations with the compound half-normal detection function.

Close modal

3. SECR approach

A total of 275 1-minute snapshots contain Cuvier's beaked whale detections on at least 1 instrument, resulting in 275 unique click-positive minutes as 275 SECR sessions (Table II). The compound half-normal model (Fig. 6) gives a much better fit than the half-normal model (ΔAIC = 47.2). In addition, estimates from the half-normal are sensitive to the choice of buffer size, likely because the data did not support the half-normal model. Hence, the compound half-normal is used for inference. With this model, the g(0) is 1.0 (jackknife CV = 0.048); the estimated detection probability remains near certain to about 1.8 km, and then drops steeply at greater distances and approaches zero by 4 km. The SECR estimate of the pseudo-EDR for a single instrument (pEDR) is 2.4 km (jackknife CV = 0.17). The mean effective survey area of the drifting array across all 275 sessions is 30.9 km2 (jackknife CV = 0.29).

TABLE II.

Summary of capture data used in the SECR analysis. Click positive minutes is the number of one-minute snapshots containing at least three Cuvier's beaked whale clicks for each recorder and the sum of all recorders. Unique click positive minutes are all one-minute snapshots with three Cuvier's beaked whale clicks on at least one DASBR and made up each session of the SECR analysis. The capture frequency is the distribution of how many unique click positive minutes were detected on either one, two, three, or all four DASBRs.

Click positive minutesUnique click positive minutesCapture frequency
B1B2B3B4ALL1234
147 195 140 162 644 275 84 73 58 60 
Click positive minutesUnique click positive minutesCapture frequency
B1B2B3B4ALL1234
147 195 140 162 644 275 84 73 58 60 
FIG. 6.

SECR estimates of the detection function of a single DASBR. Detection models are the half-normal (dashed) and the compound half-normal (solid).

FIG. 6.

SECR estimates of the detection function of a single DASBR. Detection models are the half-normal (dashed) and the compound half-normal (solid).

Close modal

1. Snapshot-based estimates

Approximately 1.22% of one-minute snapshot samples from individual instruments have Cuvier's beaked whale detections (Table III). This increases to 1.85% of the minutes if detection is considered on at least one of the four instruments (as used in the SECR estimates; Table III). The estimated effective survey area for single instruments ranges from 24.5 to 28.9 km2 for different fitting methods and detection models. The mean effective survey area for the collection of four instruments is 30.9 km2 from the SECR method. For the snapshot approaches, population density estimates range from 3.9 to 5.4 animals per 1000 km2 for three combinations of the estimation method and detection model (Table III).

TABLE III.

Estimates of the population density of Cuvier's beaked whales in the study area based on snapshot methods [Eq. (1)]. The methods for fitting the detection models include the trial-based, distance-sampling, and SECR approaches. The mean group size (1.9, CV = 0.07) and the probability of echolocation (0.209, CV = 0.04) are the same for all approaches and models. The effective survey area is for a single DASBR for the trial-based and distance-sampling approaches and for all four DASBRs for the SECR approach. The mean and CV for the effective survey area were calculated using a jackknife approach. The CV for the density is estimated by the delta method from the CVs of the quantities used to estimate it. Lower 95% confidence intervals (L95%CI) and upper 95% confidence intervals (U95%CI) are based on a log-normal distribution. The EDR for a single DASBR is a pseudo-EDR, calculated from the jackknife means of the estimated detection function parameters.

Proportion of minutes with detections n/kEDR (m) for a single DASBREffective survey area ν (km2)Density D per 1000 km2
Fitting methodDetection modelMeanCVMeanJackknife CVMeanCVL95%CIU95%CI
Trial-based Nonparametric 0.0122 0.32 2,791 24.47 0.40 4.54 0.52 1.74 11.87 
Distance-sampling Compound half-normal 0.0122 0.32 3,000 28.28 0.29 3.93 0.44 1.71 9.01 
SECR Compound half-normal 0.0185 0.37 2,376 30.90 0.29 5.44 0.48 2.23 13.28 
Proportion of minutes with detections n/kEDR (m) for a single DASBREffective survey area ν (km2)Density D per 1000 km2
Fitting methodDetection modelMeanCVMeanJackknife CVMeanCVL95%CIU95%CI
Trial-based Nonparametric 0.0122 0.32 2,791 24.47 0.40 4.54 0.52 1.74 11.87 
Distance-sampling Compound half-normal 0.0122 0.32 3,000 28.28 0.29 3.93 0.44 1.71 9.01 
SECR Compound half-normal 0.0185 0.37 2,376 30.90 0.29 5.44 0.48 2.23 13.28 

2. Dive-cue based estimates

During the 2-week study, Cuvier's beaked whale detections defined 29 distinct foraging dives (Barlow et al., 2018). The numbers of these dives as detected by the four DASBRs used in this study are 15, 21, 17, and 19. The recording effort varies between 246 and 249 h for the four instruments. The resulting overall dive detection rate is 0.073 dives hr−1 (CV = 0.24; Table IV). Assuming that the detection function for a dive is the same as the detection probability of a one-minute snapshot results in a density estimate of 15.4 animals per 1000 km2 (Table IV). Assuming that each 1-minute period within an average 39-minute dive has an independent probability of detection, the density estimate is 6.7 animals per 1000 km2 (Table IV). Ignoring the potential for movement bias, these values are expected to bracket the true dive detection probabilities, and the true density is, therefore, expected to be within this range of estimates.

TABLE IV.

Estimates of the population density of Cuvier's beaked whales in the study area based on the dive cue-counting approach [Eq. (2)] with trial-based estimates of the detection probability. The estimates include two methods of calculating the EDR (see text). The mean group size (1.9, CV = 0.07) and the mean dive rate (0.319 hr−1, CV = 0.13) are the same for both methods. The CV and mean for the effective survey area were calculated using a jackknife approach. The CV for the density is estimated by the delta method from the CVs of the quantities used to estimate it. L95%CI and U95%CI are based on a log-normal distribution.

Number of dives detected hr−1 (nd/T)Effective survey area ν (km2)Density D per 1000 km2
Fitting methodMeanCVEDR (m)MeanJackknife CVMeanCVL95%CIU95%CI
Distance-sampling, 1-minute detection probability 0.0733 0.24 3000 28.28 0.29 15.44 0.41 7.19 33.14 
Distance-sampling, 39-minute detection probability 0.0733 0.24 4539 64.73 0.16 6.74 0.32 3.64 12.51 
Number of dives detected hr−1 (nd/T)Effective survey area ν (km2)Density D per 1000 km2
Fitting methodMeanCVEDR (m)MeanJackknife CVMeanCVL95%CIU95%CI
Distance-sampling, 1-minute detection probability 0.0733 0.24 3000 28.28 0.29 15.44 0.41 7.19 33.14 
Distance-sampling, 39-minute detection probability 0.0733 0.24 4539 64.73 0.16 6.74 0.32 3.64 12.51 

By using three alternative approaches to estimating the detection functions, differences in the estimates can be compared among the approaches (Table III). For a single instrument, the EDR is 2.8 km for the trial-based approach and 3.0 km for the best-fitting compound half-normal model in the distance-sampling approach. For SECR, the equivalent pEDR for a single instrument is less (2.4 km). Although these differences in estimates of detection ranges are relatively small given their CVs, it is still informative to explore the source of these differences by examining the effect of assumptions and potential biases for each. All three estimates are based on empirical approaches that may be specific to the sound propagation conditions and animal behavior in our study area. It would not be appropriate to assume the same detection functions for other areas. These empirical approaches lead to more robust density estimates if data are used from the survey area to estimate the detection functions.

1. Trial-based

The trial-based approach is unique in that it uses the information when an animal is both heard and not heard to estimate detection probability. The primary assumptions of the trial-based approach are that the detection probability for the sample of whales at known distances is representative of the population and distances are estimated without error. This method does not assume that animals are randomly distributed in space, detection probability declines monotonically, or the probability of detection at zero horizontal range is 1.0. Also, because this method uses a smoothing spline, it does not constrain the detection function to a parametric shape. However, because our estimates are based on a small sample size (those ten animals that were tracked in the study by Barlow et al., 2018), there is a greater chance of a sampling error in applying the results to a larger set of acoustic detections.

All the tracked locations are believed to represent solitary individuals (because tracking was not feasible for larger groups; Barlow et al., 2018). Groups are more likely to be detected at a given range than a solitary individual would be because in a group, the likelihood that the relatively narrow echolocation beam of at least one individual would be pointing toward a hydrophone would be greater. If individuals moved independently when foraging as a group, the detection probability of a group of a given size could be estimated from the detection probability of an individual. Unfortunately, studies that tracked multiple individuals (Gassmann et al., 2015) show that although individuals are typically separated by several hundred meters when foraging, their overall movement is coordinated and cannot be considered independent. Detection probabilities for groups of individuals moving independently are expected to be higher than they are for solitary individuals.

Another potential predisposition in the trial-based method is that trackable individuals may be a biased subset of all individuals. Most of the tracked individuals were headed toward our group of hydrophones. By using a continuous tracking model of animal positions and multiple hydrophones, observed data include a wide distribution of animal movement relative to a potential receiving hydrophone (Fig. 2). Thus, some of this potential bias was removed by modeling the effect of orientation and averaging over all possible horizontal angles relative to a hydrophone. Although some animals were detected when their net direction of movement was away from the hydrophone, animals moving consistently away from the hydrophone array do not provide enough observations to develop a tracking model. Using a biased sample of more easily tracked individuals may result in overestimating detection probabilities (the opposite of the potential underestimation caused by including only individuals rather than groups).

Errors in estimating ranges can also bias point-transect density estimates (Buckland, 2006). Because the effective survey area is proportional to the square of the EDR, an unbiased error in the estimating range results in a biased overestimate of the area. Standard errors in tracking locations used in this study range up to 0.2 km for distant whales (see the supplementary material in Barlow et al., 2018). This leads to an overestimate of the EDR and an underestimate of the density.

The trial-based approach offers several advantages over the other two approaches in estimating a detection function. This method does not require the assumptions that animals are uniformly distributed or all animals are continuously available during the echolocation portion of their dives, and assumptions are less restrictive then for distance-sampling or SECR. However, in part due to the small sample size of tracked individuals, the CV of the effective area surveyed was higher for this method than for other methods. Only one track included observations at the greatest range (>3.8 km; Fig. 2), and the jackknife sample that excluded that track estimated a much higher EDR (3.3 km), which contributed substantially to the higher CV for this approach.

Potential biases in the trial-based method could be minimized in future studies by using a larger array of recorders. Having a wider array would provide more information on missed detections at a greater range and might thereby stabilize the estimates and reduce the associated CV. A larger array may also increase sample size and reduce the localization error and its associated bias. The potential bias caused by not being able to track multiple animals in a group is harder to address. Gassmann et al. (2015) tracked multiple independent beaked whales in a group using precisely time-aligned four-channel recording systems with a spatial hydrophone array of bottom-moored instruments. Although that configuration is difficult to replicate with a free-floating system, it may be possible to use such a system to track multiple individuals in a group and introduce the addition of drifting instruments to a bottom-moored tracking system to better quantify the detection probability of the drifting instruments as a function of the group size. If multiple individuals can be tracked in future studies, it will be important to model the detection probability as a function of the range to the center of the group to meet point-transect sampling assumptions.

2. Distance-sampling

Because distance-sampling uses only information on positive detections to estimate a detection function, more assumptions are required than for the trial-based approach. The usual assumptions are that animals are distributed randomly in horizontal space, the probability of detecting individuals at zero horizontal distance is certain (or can be estimated from other sources), and the horizontal distance is estimated without error. We modify the conventional distance-sampling approach considerably to match the available data from our acoustic study. Because the horizontal range could not be measured directly, we instead adapt this approach to use vertical detection angles measured from near-surface hydrophones. Information on horizontal distance is inferred from distributions of the detection angles (from this study) and echolocation depths (from another study). Detection probability is modeled as a function of a whale's absolute distance from a hydrophone (slant range) rather than the horizontal distance. By using a slant-range detection function, the assumption that a whale would be detected with certainty at zero range can be met, which would not necessarily be true when using the horizontal range because an animal at zero horizontal range may still be >1–2 km below the hydrophone. The simulated likelihood allows the use of maximum likelihood methods to fit a detection function in this case, where the likelihood is complex. Whereas full likelihood methods for incorporating a depth distribution in distance-sampling have been developed for the half-normal model (Thomas and Matias, 2014), the simulated likelihood is particularly flexible, allowing us to easily include the angular measurement error in our parameter estimation. Incorporating the angular error in the estimation procedure may be particularly important because the angular error has a much greater effect on the inferred range at greater distances. The distance-sampling assumptions for our modified approach are that echolocating whales can be detected with certainty at zero absolute distance (slant range), the depth distribution of echolocating whales is known without error, and the angular measurement is unbiased and its error has been accurately measured. The functional form of our parametric detection model (compound or simple half-normal) adds the further assumption that the detection function monotonically decreases and can be approximated with these functions.

Although the distance-sampling approach has more assumptions than the trial-based approach, this analysis was not limited to the subset of acoustic detections that could be tracked and, therefore, is not subject to the potential censoring biases identified for the trial-based approach. The acoustic data used to fit the distance-sampling detection function included groups of multiple individuals and a larger sample size than the trial-based approach. Perhaps most importantly, our distance-sampling approach does not require an array of instruments (which must be maintained daily) to estimate a detection function.

3. SECR

The conditional likelihood SECR approach differs from the other two approaches by not requiring an independent estimate of the distance from a hydrophone to the sound source and not requiring estimates of the vertical angle (and, consequently, not requiring two vertically spaced hydrophones). This method assumes that animals are distributed randomly in horizontal space and source locations of the calls are independent. Previous work has shown that the SECR is quite robust to nonuniform animal distribution (Efford et al., 2009a). The independence assumption did not hold true in this work as subsequent minutes containing clicks were, for some snapshots, known to be tracked individuals. However, it has been shown that violating the independence assumption leads to negligible bias in density estimates (Stevenson et al., 2015). Because nonindependence can lead to underestimation of the variance (Stevenson et al., 2015), the jackknife approach was used to estimate the variance. Additionally, acoustic SECR assumes that a detected call can be identified across hydrophones, the calls are omnidirectional, and the mean source levels do not differ by individual (Efford et al., 2009b). Although individual echolocation pulses are highly directional, we are essentially assuming that animals turn often enough during a one-minute snapshot that the detection probability is omnidirectional. Whale group densities were low in our study area, but in areas with higher densities, the assumption that each detection represented a single group would need to be evaluated.

Work has been done to improve precision and reduce bias in the parameter estimates when using SECR with acoustic data by incorporating information on signal strength, time difference of arrival, and bearing (Efford et al., 2009b; Stevenson, 2016; Stevenson et al., 2015). If signal strength information is available, directionality can also be modeled and accounted for (Stevenson, 2016). None of those improvements were used in this study because the hydrophones of each drifter were not time synchronized, and horizontal bearing information to the signals was not available. Further, because a snapshot approach was used (to overcome slow instrument movement), individual clicks (and their received levels) are not the object of interest and so signal strength information could not be incorporated in the model. Simulations demonstrate that the spacing of the four recorders and their movements during the survey are sufficient to reasonably estimate the effective survey area for these data but are likely not ideal (Fregosi, 2020). A more widely spaced array (∼1.5 km between recorders) would likely have provided increased precision of the parameter estimates and is recommended for future SECR studies of beaked whales.

Like the distance-sampling approach, the SECR analysis is not limited to the subset of acoustic detections that could be tracked and, therefore, is also not subject to the potential censoring biases identified for the trial-based approach, but it does require an array of instruments.

4. Dive detection

For the dive-based cue counting density estimates, the probability of detecting a dive as a function of the distance is needed. Our approach, assuming that the probability of detection for every minute in a typical 39-minute echolocation period is independent, is likely to overestimate the EDR if, as expected, detection is positively correlated between minutes. This 39-minute detection function is estimated only as a preliminary exploration of an alternative to the snapshot approach. The EDR for 39 independent 1-minute snapshots (4.5 km) is considerably greater than for a 1-minute snapshot (3.0 km). There are several impediments to estimating the actual detection function for a dive. First, animal movement within a 39-minute time window is not trivial. Barlow et al. (2018) estimate that the net travel speed during tracked dives used in this study was 0.63 m s−1 or approximately 1.5 km in 39 min. The dive tracks of Cuvier's beaked whales are also on the order of 1–2 km in another southern California study (Fig. 5 in Gassmann et al., 2015). Dive detection probability is likely to depend mostly on the closest distance to the hydrophone during a dive, but using that range to estimate the effective detection probability would result in a biased estimate of the population density. A simulation approach using a group movement model (Stevenson et al., 2015) may be required to accurately estimate the dive detection probability or to at least quantify the likely bias caused by movement.

5. Comparison to acoustic modeling

Most previous estimates of detection functions for beaked whales are based on acoustic models which include the echolocation pulse source level, echolocation beam patterns, propagation losses, ambient noise levels, detector characteristics, and a simulation of an animal's acoustic search behavior and movement. Many details of beaked whale searching behavior (e.g., head turning behavior as animals scan for prey with their narrow echolocation beams) are unknown and have not been included in these models, but these details could also affect the detection probability. Different authors have modeled detection probabilities for individual echolocation clicks (Küsel et al., 2011; Marques et al., 2009), time snapshots (Hildebrand et al., 2015), and entire dives (Zimmer et al., 2008). Detection probabilities for individual echolocation clicks are typically very small and not generally comparable to our estimates. Hildebrand et al. (2015) estimate that the mean probability of detecting a group of Cuvier's beaked whales within a radius of 4 km of a bottom-moored recorder in a time window of 5 min is 0.359, which corresponds to an EDR of 2.4 km. The Zimmer et al. (2008) model predicts that the dive detection probability for a solitary Cuvier's beaked whale using a single hydrophone at 100 m depth depends on whether the animal periodically reversed direction during the course of a dive. If it continues straight during a dive, the detection probability was ∼1.0 until 0.7 km, dropped to ∼0.5 at 2 km, and then dropped to zero at 4 km. If animals turn during a dive, the detection probability remained near 1.0 to a range of 2 km and then dropped steeply to zero at 4 km. Zimmer et al. (2008) does not estimate the effective detection radii for their models, but based on the similarity in shapes, their model with periodic direction reversals would have an EDR for dive detections that would be very similar to those of our compound half-normal model fit by distance-sampling (3.0 km). These acoustic models do not consider all sources of uncertainty and most are for seafloor hydrophones, therefore, it is not possible to statistically compare our empirically based estimates of the detection probability; however, our empirical estimates of the detection range are near the upper ranges of estimates based on the acoustic modeling.

6. Detection probability at zero distance: g(0)

All three empirical approaches to estimating detection functions estimate a nearly certain probability of detection at zero horizontal distance. However, the data used in this study had no observations that were closer than 0.5 km to a hydrophone (Fig. 2), hence, this result should be considered with caution. An acoustic “shadow” could exist with a lower probability of detection for animals that are directly below a hydrophone. This would be expected if animals primarily search horizontally with their echolocation “beams” and seldom search straight up. There is little data on the details of beaked whale foraging behavior, but Laplanche et al. (2015) modeled the 3D diving behavior of a tagged Blainville's beaked whale and found that its pitch angle during foraging was centered at 0° (relative to horizontal) and seldom exceeded ±50° (Laplanche et al., 2015, supplemental appendix S5). Our trial-based study shows a high probability of detection at ranges of 0.5–1.0 km (Fig. 2). In point-transect sampling, a lowered probability of detection near the sampling location has a small effect on the EDR because the sampled area close to that point is small compared to the area at greater distances. Nonetheless, it would be good to have more data at close range to evaluate whether a shadow effect may be real.

7. Recommendations

Because distance-sampling can estimate a detection function with data from a single instrument, this is clearly the most efficient and practical approach for future large-scale acoustic surveys of beaked whale density and abundance. The SECR approach requires data from an array of at least four instruments, and the trial-based approach requires data from a tracking model that was developed using an array of five instruments. For a short study in a small study area, it was practical to reestablish the array configuration needed for the SECR and tracking on a daily basis. This would not be the case with an acoustic survey over greater temporal and spatial scales.

Although simultaneous use of three approaches to estimate the EDR may be impractical for large-scale acoustic surveys, the use of all three in one area is informative. The distance-sampling approach depends critically on some key assumptions which are, to some extent, validated by the other two approaches. All methods estimated a very high probability of detection at ranges up to 1.5 km. The EDR estimates from the distance-sampling approach with the compound half-normal model (3.0 km) are generally consistent with the estimate from the trial-based approach (2.8 km) given that the latter is likely to be biased low by only including solitary individuals. It is not clear why the pEDR from the SECR (2.4 km) is less than the EDR of the other two; however, the measurement error could explain at least part of the observed differences among all the approaches. A better understanding of the differences among the methods could be obtained with a larger sample size collected from a larger array of instruments.

1. Availability bias

Snapshot-based point-transect estimates of density are sensitive to availability bias when using our empirically estimated detection functions. This is addressed in our density formula [Eq. (1)] by including a term that corrects for times when animals are not foraging and producing echolocation pulses. However, whales may not be detectable during the entirety of the echolocation periods, and periods of unavailability could result in an underestimate of abundance if they last longer than the time period of a snapshot. In principle, our trial-based method should correct for short periods of unavailability because that method uses information on missed detections; however, our sample of tracked individuals may be a biased subset of dives that did not include long periods of unavailability. Intermittent availability during the echolocation portion of a dive could result from a variety of factors. During the initial descent phase of a foraging dive, beaked whales descend steeply from the surface until they reach their foraging depth (Tyack et al., 2006; Laplanche et al., 2015, supplemental appendix S5). As an animal descends vertically from when the echolocation starts (∼500 m) until it reaches its foraging depth (>800 m), its highly directional echolocation beam is not likely to be detected at a near-surface hydrophone unless the absolute distance to the hydrophone is less than 700 m (Zimmer et al., 2008). Also, during foraging, an animal might be intermittently unavailable for detection when it is headed away from a hydrophone. The shadow effect (low detectability for animals directly below the hydrophone) may contribute to intermittent availability. Overall, the fraction of time available for detection is likely overestimated by assuming availability during the entirety of their echolocation period and, thus, the beaked whale density may be underestimated. Because the distance-sampling and SECR approaches are based only on detections, their estimated detection functions only apply to times when animals are available to be detected. The trial-based approach is less susceptible to this type of bias because it uses data both from times when animals are detected and when they are not. Trial-based detection functions, thus, include an element of availability, which might also explain why trial-based estimates of the EDR are lower than for the other two approaches.

Availability bias in the snapshot approach can be reduced by using longer snapshots (increasing the likelihood that a group will be available at some point during the snapshot); however, animal movement within a longer time snapshot can also bias abundance estimates (Buckland et al., 2015). Additional work is needed to determine the optimum time window for use with the snapshot approach to minimize the potential biases from movement and intermittent availability.

In estimating availability, we have assumed that individual whales in a group start and end echolocation synchronously during a dive. Aguilar de Soto (2018) found beaked whale group dives to be highly synchronous with a 98% overlap in the echolocation period when two animals in the same group were tagged with acoustic recording tags. Asynchrony would increase the availability by increasing the time available for group detection, but this effect appears to be small.

2. Snapshot vs dive cue methods

Both the snapshot and dive cue methods of density estimation depend on information that is external to an acoustic survey. Both depend on estimates of the group size, which are taken here from visual surveys. The snapshot density estimates also depend on the fraction of time in which whales are available to be detected, which is estimated from tagging studies as the fraction of time at foraging depths (Barlow et al., 2020). The dive cue-count density estimates also depend on the dive rate estimated from the tagging studies (Barlow et al., 2020). In this study, all these external sources of information were collected in the vicinity of our study area. This may not be the case for other studies, and there is a need to quantify the variability of these parameters between study sites if our methods are to be applied more widely (Warren et al., 2017). Sampling-based density estimates are most robust if all the needed parameters can be estimated empirically within a single survey (Buckland et al., 2015). Although it may never be practical to estimate availability (the proportion of time echolocating) from an acoustic-only survey, it may be possible to empirically estimate the depth distribution of echolocating whales from a single instrument using surface reflections (Barlow and Griffiths, 2017). In some cases, it may also be possible to estimate dive rates from the time period between the echolocation signals received from successive dives by the same group.

The dive cue-counting approach has an advantage over the snapshot-based approach because the assumption that all dives are detected at zero horizontal distance should be easier to meet than the assumption that all one-minute snapshots are detected. Intermittent availability would also be less of a problem for dive cue-counting because there are many more opportunities to detect a group during a dive than there are during a one-minute snapshot. Our dive-based density estimate is higher than any of the snapshot-based estimates, which is consistent with expectations if the snapshot estimates are biased low by intermittent availability. The CV of the dive-based density estimate is also lower than any of the other approaches. At this time, given the uncertainty in estimating an acoustic detection function for dives and the potential for bias caused by animal movement during a dive (Glennie et al., 2015), the dive cue-count method cannot be considered reliable. A full consideration of variation in dive lengths is also needed (instead of simply using a mean dive duration). Additional research is needed, but this approach has great potential.

In contrast, there are few impediments to the immediate use of the snapshot-based method to estimate beaked whale densities. Using the distance-sampling approach, the primary sources of information for this method (depth distributions of echolocation, fraction of time echolocating, and mean group size) are well quantified for Cuvier's beaked whales. Less external information (outside of what can be collected on a survey) is needed for the empirical approaches developed here than for previous acoustic modeling approaches. Because estimates of dive rates or acoustic availability are needed to apply these methods to other beaked whale species, tagging studies should be expanded to include more species.

3. Comparison to previous density estimates

Our acoustic-based estimates of Cuvier's beaked whale density (3.9–6.7 animals per 1000 km2) are higher than a previous visual survey estimate for the entire U.S. west coast study area (3.2 animals per 1000 km,2 CV = 0.35; calculated from Barlow, 2016, Tables 2 and 7). There are insufficient sighting data to make visual estimates of the beaked whale density within the tiny area surveyed in this study. Indeed, there were only 63 sightings of Cuvier's beaked whales on 7 visual sighting surveys from 1991 to 2014. Those surveys covered transect lines totaling 72 123 km, and each survey required several months on a large research ship. Our ability to detect 29 distinct dives of Cuvier's beaked whales in a 2-week period illustrates the potential power of acoustic-based survey methods to improve density and abundance estimations for this species.

The only previous estimates of Cuvier's beaked whale density from an acoustic survey in Southern California were made by Hildebrand et al. (2016). Their snapshot-based and cue-count-based mean estimates for a single location (their site N) range from 1.7 to 1.9 animals per 1000 km2, respectively. For both estimates, detection functions were estimated using an acoustic modeling approach.

Although the CVs of our density estimates are moderately high (0.32–0.52), estimates from visual sighting surveys are typically much higher. The CVs of the visual estimates of Cuvier's beaked whale abundance along the U.S. west coast range from 0.51 to 0.68 for five visual surveys from 1996 to 2014 (Table 8 in Barlow, 2016). Each of those surveys required 3–5 months of ship time and a much greater operational cost. Obtaining CVs in the range of 0.32–0.52 with two weeks of field effort is remarkable and hints at the potential benefit of passive acoustic surveys for beaked whales.

Funding for this research was provided by the U.S. Office of Naval Research (Project Nos. N00014-15-1-2142 and MIPR N0001416IP00059) and NOAA's Southwest Fisheries Science Center. Vessel time on the Horizon was funded by NOAA's Cooperative Research Program. Funding for DASBR development and some equipment was provided by the U.S. Navy's N45 and Living Marine Resource programs and NOAA's Acoustics Program. S.F. was supported by the National Science and Engineering Graduate Fellowship. We thank Mike Weise, Frank Stone, Jason Gedamke, Lisa Ballance, and Anu Kumar for their support. Field assistance was provided by Dave Mellinger, Holger Klinck, Jennifer McCullough, and Eiren Jacobson. Vessel operators were Trevor Oudin, Juan Carlos Aguilar, and Spencer Salmon. Analysis methods were developed with helpful ideas from Tiago Marques and Jeff Moore. This manuscript was improved by helpful comments from Jeff Moore. This is Pacific Marine Environmental Laboratory (PMEL) Contribution No. 5071.

1

See supplementary material at https://www.scitation.org/doi/suppl/10.1121/10.0002881 for the R code for the trial-based, distance-sampling, and SECR estimates of the detection probability, additional methods for the simulated likelihood method used in the distance-sampling approach, and model fit figures for the distance-sampling approach.

1.
Aguilar de Soto
,
N.
,
Visser
,
F.
,
Madsen
,
P. T.
,
Tyack
,
P.
,
Ruxton
,
G.
,
Alcazar
,
J.
,
Arranz
,
P.
, and
Johnson
,
M.
(
2018
). “Beaked and killer whales show how collective prey behaviour foils acoustic predators,” BioRxiv:
303743
.
2.
Ainslie
,
M. A.
(
2013
). “
Neglect of bandwidth of Odontocetes echo location clicks biases propagation loss and single hydrophone population estimates
,”
J. Acoust. Soc. Am.
134
,
3506
3512
.
3.
Amante
,
C.
, and
Eakins
,
B. W.
(
2009
). “
ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources, and Analysis
,” NOAA Technical Memorandum NESDIS NGDC-24.
4.
Arias
,
C.
, and
Cox
,
T. L.
(
1999
). “
Maximum simulated likelihood: A brief introduction for practitioners
,” in
Agricultural and Applied Economics Staff Papers Series
, Staff Paper No. 421, University of Wisconsin–Madison.
5.
Au
,
W. W. L.
(
1993
).
The Sonar of Dolphins
(
Springer
,
New York
).
6.
Barlow
,
J.
(
2016
). “
Cetacean abundance in the California Current estimated from ship-based line-transect surveys in 1991–2014
,” NOAA Southwest Fisheries Science Center Administrative Report LJ-16-01,
63
pp.
7.
Barlow
,
J.
, and
Gisiner
,
R.
(
2006
). “
Mitigating, monitoring and assessing the effects of anthropogenic sound on beaked whales
,”
J. Cetacean Res. Manag.
7
,
239
249
.
8.
Barlow
,
J.
, and
Griffiths
,
E. T.
(
2017
). “
Precision and bias in estimating detection distances for beaked whale echolocation clicks using a two-element vertical hydrophone array
,”
J. Acoust. Soc. Am.
141
,
4388
4397
.
9.
Barlow
,
J.
,
Griffiths
,
E. T.
,
Klinck
,
H.
, and
Harris
,
D. V.
(
2018
). “
Diving behavior of Cuvier's beaked whales inferred from three-dimensional acoustic localization and tracking using a nested array of drifting hydrophone recorders
,”
J. Acoust. Soc. Am.
144
2030
2041
.
10.
Barlow
,
J.
,
Schorr
,
G. S.
,
Falcone
,
E. A.
, and
Moretti
,
D.
(
2020
). “
Variation in dive behavior of Cuvier's beaked whales with seafloor depth, time-of-day, and lunar illumination
,”
Mar. Ecol. Progr. Ser.
644
,
199
214
.
11.
Barlow
,
J.
,
Tyack
,
P. L.
,
Johnson
,
M. P.
,
Baird
,
R. W.
,
Schorr
,
G. S.
,
Andrews
,
R. D.
, and
Aguilar de Soto
,
N.
(
2013
). “
Trackline and point detection probabilities for acoustic surveys of Cuvier's and Blainville's beaked whales
,”
J. Acoust. Soc. Am.
134
,
2486
2496
.
12.
Baumann-Pickering
,
S.
,
McDonald
,
M. A.
,
Simonis
,
A. E.
,
Solsona Berga
,
A.
,
Merkens
,
K. P. B.
,
Oleson
,
E. M.
,
Roch
,
M. A.
,
Wiggins
,
S. M.
,
Rankin
,
S.
,
Yack
,
T. M.
, and
Hildebrand
,
J. A.
(
2013
). “
Species-specific beaked whale echolocation signals
,”
J. Acoust. Soc. Am.
134
,
2293
2301
.
13.
Borchers
,
D.
(
2012
). “
A non-technical overview of spatially explicit capture-recapture models
,”
J. Ornithol.
152
,
435
444
.
14.
Borchers
,
D. L.
, and
Efford
,
M. G.
(
2008
). “
Spatially explicit maximum likelihood methods for capture-recapture studies
,”
Biometrics
64
,
377
385
.
15.
Buckland
,
S. T.
(
2006
). “
Point-transect surveys for songbirds: Robust methodologies
,”
Auk
123
,
345
357
.
16.
Buckland
,
S. T.
,
Rexstad
,
E. A.
,
Marques
,
T. A.
, and
Oedekoven
,
C. S.
(
2015
).
Distance Sampling: Methods and Applications
(
Springer Science and Business Media
,
New York)
.
17.
Dawson
,
S.
,
Wade
,
P.
,
Slooten
,
E.
, and
Barlow
,
J.
(
2008
). “
Design and field methods for sighting surveys of cetaceans in coastal and riverine habitats
,”
Mamm. Rev.
38
,
19
49
.
18.
Efford
,
M. G.
(
2019a
). “
Multi-session models in secr 4.1
,” available at https://www.otago.ac.nz/density/ (Last viewed 18 December 2020).
19.
Efford
,
M. G.
(
2019b
). “
secr—Spacially explicit capture recapture in R
,” available at https://www.otago.ac.nz/density/ (Last viewed 18 December 2020).
20.
Efford
,
M. G.
,
Borchers
,
D. L.
, and
Byrom
,
A. E.
(
2009a
). “
Density estimation by spatially explicit capture-recapture: Likelihood-based methods
,” in
Modeling Demographic: Processes in Marked Populations
, edited by
D. L.
Thomson
,
E. G.
Cooch
, and
M. J.
Conroy
(
Springer
,
Boston)
, pp.
255
269
.
21.
Efford
,
M. G.
, and
Dawson
,
D. K.
(
2009
). “
Effect of distance-related heterogeneity on population size estimates from point counts
,”
Auk
126
,
100
111
.
22.
Efford
,
M. G.
,
Dawson
,
D. K.
, and
Borchers
,
D. L.
(
2009b
). “
Population density estimation from locations of individuals on a passive detector array
,”
Ecology
90
,
2676
2682
.
23.
Efron
,
B.
(
1982
). “
The jackknife, the bootstrap and other resampling plans
,” in
CBMS Regional Conference Series in Applied Mathematics 38
(
Society for Industrial and Applied Mathematics
,
Philadelphia, PA
).
24.
Fregosi
,
S.
(
2020
). “
Applications of slow-moving autonomous platforms for passive acoustic monitoring and density estimation of marine mammals
,” Ph.D. thesis,
Oregon State University
,
Corvallis, OR
.
25.
Gassmann
,
M.
,
Wiggins
,
S. M.
, and
Hildebrand
,
J. A.
(
2015
). “
Three-dimensional tracking of Cuvier's beaked whales' echolocation sounds using nested hydrophone arrays
,”
J. Acoust. Soc. Am.
138
,
2483
2494
.
26.
Gillespie
,
D.
,
Mellinger
,
D. K.
,
Gordon
,
J.
,
McLaren
,
D.
,
Redmond
,
P.
,
McHugh
,
R.
,
Trinder
,
P.
,
Deng
,
X.-Y.
, and
Thode
,
A.
(
2009
). “
PAMGUARD: Semiautomated, open source software for real-time acoustic detection and localization of cetaceans
,”
J. Acoust. Soc. Am.
125
,
2547
.
27.
Glennie
,
R.
,
Buckland
,
S. T.
, and
Thomas
,
L.
(
2015
). “
The effect of animal movement on line transect estimates of abundance
,”
PloS One
10
,
e0121333
.
28.
Heinemann
,
D.
,
Gedamke
,
J.
,
Oleson
,
E.
,
Barlow
,
J.
,
Crance
,
J.
,
Holt
,
M.
,
Soldevilla
,
M.
, and
Van Parijs
,
S.
(
2016
). “
Report of the Joint Marine Mammal Commission—National Marine Fisheries Service Passive Acoustic Surveying Workshop
,” NOAA Technical Memorandum NMFS-F/SPO-164,
107
pp.
29.
Hildebrand
,
J. A.
,
Baumann-Pickering
,
S.
,
Frasier
,
K. E.
,
Trickey
,
J. S.
,
Merkens
,
K. P.
,
Wiggins
,
S. M.
,
McDonald
,
M. A.
,
Garrison
,
L. P.
,
Harris
,
D.
,
Marques
,
T. A.
, and
Thomas
,
L.
(
2015
). “
Passive acoustic monitoring of beaked whale densities in the Gulf of Mexico
,”
Sci. Rep.
5
,
16343
.
30.
Hildebrand
,
J. A.
,
Baumann-Pickering
,
S.
,
Giddings
,
A.
,
Brewer
,
A.
,
Jacobs
,
E. R.
,
Trickey
,
J. S.
,
Wiggins
,
S. M.
, and
McDonald
,
M. A.
(
2016
). “
Progress report on the application of passive acoustic monitoring to density estimation of Cuvier's beaked whales
,” MPL Technical Memorandum No. 606,
16
pp.
31.
Keating
,
J. L.
, and
Barlow
,
J.
(
2013
). “
Click detectors and classifiers used during the 2012 Southern California Behavioral Response Study
,” NOAA Technical Memorandum NOAA-TM-NMFS-SWFSC-517,
14
pp.
32.
Küsel
,
E. T.
,
Mellinger
,
D. K.
,
Thomas
,
L.
,
Marques
,
T. A.
,
Moretti
,
D.
, and
Ward
,
J.
(
2011
). “
Cetacean population density estimation from single fixed sensors using passive acoustics
,”
J. Acoust. Soc. Am.
129
,
3610
3622
.
33.
Laplanche
,
C.
,
Marques
,
T. A.
, and
Thomas
,
L.
(
2015
). “
Tracking marine mammals in 3D using electronic tag data
,”
Methods Ecol. Evol.
6
,
987
996
.
34.
Marques
,
T. A.
,
Jorge
,
P. A.
,
Mouriño
,
H.
,
Thomas
,
L.
,
Moretti
,
D. J.
,
Dolan
,
K.
,
Claridge
,
D.
, and
Dunn
,
C.
(
2019
). “
Estimating group size from acoustic footprint to improve Blainville's beaked whale abundance estimation
,”
Appl. Acoust.
156
,
434
439
.
35.
Marques
,
T. A.
,
Munger
,
L.
,
Thomas
,
L.
,
Wiggins
,
S.
, and
Hildebrand
,
J. A.
(
2011
). “
Estimating North Pacific right whale Eubalaena japonica density using passive acoustic cue counting
,”
Endanger. Species Res.
13
(
3
),
163
172
.
36.
Marques
,
T. A.
,
Thomas
,
L.
,
Martin
,
S. W.
,
Mellinger
,
D. K.
,
Jarvis
,
S.
,
Morrissey
,
R. P.
,
Ciminello
,
C.-A.
, and
DiMarzio
,
N.
(
2012
). “
Spatially explicit capture-recapture methods to estimate minke whale density from data collected at bottom-mounted hydrophones
,”
J. Ornithol.
152
,
445
455
.
37.
Marques
,
T. A.
,
Thomas
,
L.
,
Martin
,
S. W.
,
Mellinger
,
D. K.
,
Ward
,
J. A.
,
Moretti
,
D. J.
,
Harris
,
D.
, and
Tyack
,
P.
(
2013
). “
Estimating animal population density using passive acoustics
,”
Biol. Rev.
88
,
287
309
.
38.
Marques
,
T. A.
,
Thomas
,
L.
,
Ward
,
J.
,
DiMarzio
,
N.
, and
Tyack
,
P. L.
(
2009
). “
Estimating cetacean population density using fixed passive acoustic sensors: An example with Blainville's beaked whales
,”
J. Acoust. Soc. Am.
125
,
1982
1994
.
39.
Mathias
,
D.
,
Thode
,
A. M.
,
Straley
,
J.
, and
Andrews
,
R. D.
(
2013
). “
Acoustic tracking of sperm whales in the Gulf of Alaska using a two-element vertical array and tags
,”
J. Acoust. Soc. Am.
134
,
2446
2461
.
40.
McDonald
,
M. A.
, and
Fox
,
C. G.
(
1999
). “
Passive acoustic methods applied to fin whale population density estimation
,”
J. Acoust. Soc. Am.
105
,
2643
2651
.
41.
Moretti
,
D.
,
Marques
,
T. A.
,
Thomas
,
L.
,
DiMarzio
,
N.
,
Dilley
,
A.
,
Morrissey
,
R.
,
McCarthy
,
E.
,
Ward
,
J.
, and
Jarvis
,
S.
(
2010
). “
A dive counting density estimation method for Blainville's beaked whale (Mesoplodon densirostris) using a bottom-mounted hydrophone field as applied to a mid-frequency active (MFA) sonar operation
,”
Appl. Acoust.
71
,
1036
1042
.
42.
R Core Team
(
2017
). “
R: A language and environment for statistical computing
,” R Foundation for Statistical Computing, Vienna, Austria, available at https://www.R-project.org/ (Last viewed 18 December 2020).
43.
Royle
,
J. A.
,
Chandler
,
R. B.
,
Sollmann
,
R.
, and
Gardner
,
B.
(
2013
).
Spatial Capture-Recapture
(
Elsevier
,
Boston, MA)
.
44.
Schorr
,
G. S.
,
Falcone
,
E. A.
,
Moretti
,
D. J.
, and
Andrews
,
R. D.
(
2014
). “
First long-term behavioral records from Cuvier's beaked whales (Ziphius cavirostris) reveal record-breaking dives
,”
PLoS One
9
,
e92633
.
45.
Stevenson
,
B. C.
(
2016
). “
Methods in spatially explicit capture-recapture
,” Ph.D. thesis,
University of St. Andrews
,
St. Andrews, Scotland
.
46.
Stevenson
,
B. C.
,
Borchers
,
D. L.
,
Altwegg
,
R.
,
Swift
,
R. J.
,
Gillespie
,
D. M.
, and
Measey
,
G. J.
(
2015
). “
A general framework for animal density estimation from acoustic detections across a fixed microphone array
,”
Methods Ecol. Evol.
6
,
38
48
.
47.
Thomas
,
L.
, and
Matias
,
L.
(
2014
). “
Cheap DECAF: Density Estimation for Cetaceans from Acoustic Fixed sensors using separate, non-linked devices
,” Final report for the Office of Naval Research. Available from contact information posted on: https://www.onr.navy.mil/Science-Technology/Departments/Code-32/all-programs/marine-mammals-biology (Last viewed 21 December 2020).
48.
Tyack
,
P. L.
,
Johnson
,
M.
,
Aguilar de Soto
,
N.
,
Sturlese
,
A.
, and
Madsen
,
P. T.
(
2006
). “
Extreme diving of beaked whales
,”
J. Exp. Biol.
209
,
4238
4253
.
49.
Urian
,
K.
,
Gorgone
,
A.
,
Read
,
A.
,
Balmer
,
B.
,
Wells
,
R. S.
,
Berggren
,
P.
,
Durban
,
J.
,
Eguchi
,
T.
,
Rayment
,
W.
, and
Hammond
,
P. S.
(
2015
). “
Recommendations for photo-identification methods used in capture-recapture models with cetaceans
,”
Mar. Mammal Sci.
31
,
298
321
.
50.
Warren
,
V. E.
,
Marques
,
T. A.
,
Harris
,
D.
,
Thomas
,
L.
,
Tyack
,
P. L.
,
Aguilar de Soto
,
N.
,
Hickmott
,
L. S.
, and
Johnson
,
M. P.
(
2017
). “
Spatio-temporal variation in click production rates of beaked whales: Implications for passive acoustic density estimation
,”
J. Acoust. Soc. Am.
141
,
1962
1974
.
51.
Wiggins
,
S. M.
,
McDonald
,
M. A.
,
Munger
,
L. M.
,
Moore
,
S. E.
, and
Hildebrand
,
J. A.
(
2004
). “
Waveguide propagation allows range estimates for North Pacific right whales in the Bering Sea
,”
Can. Acoust.
32
,
146
154
.
52.
Wood
,
S. N.
(
2017
).
Generalized Additive Models: An Introduction with R
(
CRC Press
,
Boca Raton, FL
).
53.
Wood
,
S.
, and
Scheipl
,
F.
(
2017
). “
gamm4: Generalized Additive Mixed Models using ‘mgcv’ and ‘lme4
,’” R package version 0.2-5, available at https://CRAN.R-project.org/package=amm4 (Last viewed 18 December 2020).
54.
Zimmer
,
W. M. X.
,
Harwood
,
J.
,
Tyack
,
P. L.
,
Johnson
,
M. P.
, and
Madsen
,
P. T.
(
2008
). “
Passive acoustic detection of deep-diving beaked whales
,”
J. Acoust. Soc. Am.
124
,
2823
2832
.

Supplementary Material