Abstract
Temporal variability is an important feature of aquatic ecosystems that can be difficult to measure. A stationary, upward-facing scientific echosounder was used to record the vertical distribution of pelagic fauna in Monterey Bay, CA, for 18 months. To characterize the distributions, a suite of metrics, including measures of density, abundance, location, dispersion, occupancy, evenness, and aggregation, was developed and tested. An algorithm to detect and count the number of acoustic backscatter layers was developed using image-analysis techniques. The metrics recorded a strong seasonal cycle, with total backscatter reaching a minimum during the spring upwelling season and peaking in autumn and winter. Variability in the vertical distribution of animals was greatest at long time-scales and decreased as a power (−1.050 to −1.585) of signal frequency. There were significant peaks in the power spectrum at 12- and 24-h periods, corresponding to the semi-diurnal tide and diel vertical migration. The diel signal was strongest in late winter and weakest during the spring upwelling season. Active acoustics are a useful addition to ocean observatories, and the metrics presented provide a useful set of tools to quantify the distribution and temporal variability of pelagic fauna.Urmy, S. S., Horne, J. K., and Barbee, D. H. 2012. Measuring the vertical distributional variability of pelagic fauna in Monterey Bay. – ICES Journal of Marine Science, 69: 184–196.
Introduction
Aquatic ecosystems display dynamic behaviour across a wide range of temporal and spatial scales, with physical and biological processes occurring both quickly and slowly over short and long distances (Stommel, 1963; Haury et al., 1978). Understanding these processes depends on understanding their variability. Temporal variability has been measured less frequently than spatial variability, especially in aquatic ecosystems. Long-term, high-resolution measurements are common for physical variables, but biologists have typically had to be content with time-series that are either short with high resolution or long with a low sampling frequency. In an incompletely understood system, the scope of observations (i.e. their extent over the smallest sample size) determines their power to detect patterns present at multiple scales (Schneider, 2009). If true, then the relative scarcity of high-scope biological time-series may represent a significant gap in our ecological knowledge.
The most common approach used to describe and quantify marine environments has been ship-based, mobile surveys. Constraints of the ship-based approach include the difficulty of collecting high-scope time-series and the confounding of space and time when traversing a study area. An alternative approach is marine observatories: typically stationary platforms, on the surface or bottom, providing an attachment point, power, and/or a communications link for instruments. Ocean observatories complement ship-based research and are poised to play an increasing role in marine biology. All ocean observatories measure physical variables, at the very least, temperature and salinity. Many measure chlorophyll or fluorescence as proxies for primary production, but few are equipped to monitor higher trophic levels. Active acoustics (i.e. sonar) is a technology capable of remotely sensing the distribution and density of animals underwater. Traditionally deployed on ships in conjunction with net-sampling, scientific echosounders are also suited for observatory use, providing a synoptic view of fish and zooplankton through the water column with high resolution in both space and time. Taxonomic resolution is coarse, especially in long-term deployments, because of difficulties in acoustic species identification (Horne, 2000). Used as a proxy for biomass, analogous to satellite ocean colour for primary productivity, acoustic backscatter can be used to characterize patterns and dynamics in the pelagic zone through time.
To date, there have been relatively few deployments of bottom-mounted and moored echosounders, with most studies focused on the behaviour of fish in freshwater (e.g. Arrhenius et al., 2000; Prchalová et al., 2003; Mehner, 2006) and the coastal ocean (e.g. Axenrot et al., 2004; Kaartvedt et al., 2009; Kaltenberg and Benoit-Bird, 2009). Farther offshore, stationary acoustics have been used to monitor micronekton (Trevorrow, 2005) and marine mammals (Doksæter et al., 2009). Others have used stationary acoustics to investigate specific ecological processes, including bird foraging on mesopelagic animals in pack ice (Kaufmann et al., 1995), copepod diapause (Osgood and Checkley, 1997), and the effect of physical processes on thin layers of zooplankton (McManus et al., 2005).
A few studies have taken a more general approach, using active acoustics to characterize the aggregate state and variability of the pelagic ecosystem over multiple months and years. Flagg et al. (1994) used acoustic Doppler current profilers (ADCPs) to record currents and biological scattering, in conjunction with chlorophyll a measurements over 15 months in the Mid-Atlantic Bight. Cochrane et al. (1994) used an ADCP to measure the temporal variability of euphausiid (Meganyctiphanes norvegica) concentrations over 49 d on the Scotian shelf. Brierley et al. (2006) used an acoustic mooring to monitor changes in the abundance of Antarctic krill (Euphausia superba) over 3 months at South Georgia Island in the Southern Ocean. Radenac et al. (2010) analysed ADCP data from the TAO/TRITON array in the equatorial Pacific, linking patterns in the acoustic backscatter with environmental forcing. ADCPs are widely deployed on oceanographic moorings, but quantitative interpretation of their backscatter data is limited by difficulties with their calibration (Brierley et al., 1998).
This study used a bottom-mounted, upward-facing scientific echosounder to collect water column profiles over 18 months at the Monterey Accelerated Research System (MARS, http://www.mbari.org/mars/) observatory in Monterey Bay, CA. The objectives were twofold. The first was to develop methods to characterize the distribution of animals in the water column. To this end, a suite of summary metrics was developed to measure characteristics and features of the vertical density distribution. The second objective was to use these metrics to quantify temporal variability in Monterey Bay through time and across temporal scales, identifying dominant modes of variability and significant changes in the ecosystem.
Methods
Instruments and data
Monterey Bay is a large open embayment on the central California coast. It is mainly shallow (<100 m), but is bisected by the Monterey Submarine Canyon, which brings deep-water habitat close to shore in the centre of the Bay. The Bay is situated in the California Current System, and its oceanographic seasons are defined by the presence or the absence of wind-driven upwelling. Upwelling events are episodic during spring and summer, tailing off into a late-summer and early autumn “oceanic period”, followed by the winter appearance of the north-flowing Davidson Current (Huyer, 1983; Pennington and Chavez, 2000).
The MARS observatory node, operated by the Monterey Bay Aquarium Research Institute (MBARI), is located at 36°42.75′N 122°11.21′W near 900-m depth on the continental slope north of the Canyon (Figure 1). MARS is connected to shore by a 52-km cable, providing continuous power and communications for up to eight scientific instruments. The node consists of a 3.7 × 4.6 × 1.2 m trawl-resistant metal frame housing an electronics assembly, which regulates power and transmits data to shore. Instruments are deployed and recovered by remotely operated vehicles (ROVs) from MBARI vessels.
Monterey Bay, showing the MARS node (black square) and the cable (black line). Isobaths are at 500 m intervals. The inset shows the location on the coast of California. Bathymetry after Carnigan et al. (2009).
The Deep Echo Integrating Marine Observatory System (DEIMOS) was deployed at −875 m from 27 February 2009 to 18 August 2010. An active acoustic observing package, DEIMOS, was built around a Simrad EK60 38 kHz scientific echosounder. The transducer had a beam width of 7° (between half-power points) and was oil-filled, allowing operation to depths of 1500 m. The transceiver and other electronics were housed in a borosilicate glass pressure sphere and mounted, with the transducer, on a galvanized steel frame. The echosounder sampled the water column at a frequency of 0.2 Hz, using a pulse length of 1.024 ms. Transmitted power was restricted to 825 W to conform to sound pressure levels specified in the US Marine Mammal Protection Act and the Endangered Species Act.
Over the 18-month deployment, data collection was interrupted several times, because of problems with software, communications, and the power supply. The longest of these interruptions was from 18 May to 14 August 2009, when increased electrical noise saturated the DEIMOS receiver and rendered the data unusable. There were three other multiday outages: 3–13 September 2009 (software crash), 23 March–7 April 2010 (power supply noise), and 22 July–4 August 2010 (shore cable severed by burrowing rodents).
DEIMOS was calibrated in situ using a 38.1-mm tungsten-carbide reference sphere following procedures outlined in Foote et al. (1987). On 24 June 2010, the ROV “Ventana” put the calibration apparatus in place. The sphere was held 12 m above the transducer face, anchored by two small lead weights on either side of the transducer, and suspended from a syntactic foam float 5 m above. All calibration components were connected by a monofilament line. The sphere moved through the acoustic beam with the currents and completed beam coverage in 3–4 d. The calibration sphere was left in place for 8 weeks and removed on 17 August, 1 d before DEIMOS was recovered. In that period, we conducted nine separate calibrations. For each calibration, beam-angle corrections were calculated using Simrad's LOBE program, and transducer gains and area-backscattering coefficient (Sa) corrections were calculated by integrating all on-axis target detections using Echoview software (v. 4.9, Myriax Pty Ltd, 2010). Before analysis was started, we corrected all raw acoustic data using the mean values of these calibration parameters.
Acoustic datafiles were processed using Echoview. Initial processing removed noise interference (e.g. electrical fluctuations, passing ships) and identified data-gaps. The top 10 m and bottom 7 m of the water column were excluded from the analysis to avoid integrating bubbles from breaking waves or echoes in the acoustic nearfield. Background noise was estimated and subtracted using methods reported in DeRobertis and Higginbottom (2007), who recommended collecting data at least 600 m below the bottom (in our case, above the surface) for unbiased noise estimates. Although we only collected data to ∼300 m above the surface, we did not find any difference between the estimates from their algorithm and standard estimates using passive listening data (Nunnallee, 1990). A signal-to-noise ratio (SNR) threshold was not applied to the final de-noised data, because during extended periods of slightly higher background noise, mean volumebackscattering values were biased low by eliminating pixels with a low SNR. After elimination of bad data regions and noise subtraction, acoustic data were exported in 5 s horizontal by 0.5 m vertical matrices of mean volumebackscattering strength (Sv; MacLennan et al., 2002) for further analysis.
Metric selection and testing
To describe this large dataset parsimoniously, we developed a suite of metrics to characterize vertical distributions of acoustic backscatter through time. Acoustic backscatter is interpreted as a proxy for the density of aquatic organisms (Foote, 1983). Metrics were derived from backscatter, as well as indices used to describe the spatial structure and variance of animal population densities (Bez and Rivoirard, 2001; Woillez et al., 2007; Burgos and Horne, 2008). These metrics quantify the vertical distribution of pelagic animals, including total abundance, mean density, occupied area, mean location, spread, evenness, aggregation, and the number of backscattering layers (Table 1). In the descriptions below, volume and area-backscattering strengths (Sv and Sa) are the logarithmic forms of the volume and area-backscattering coefficients sv and sa [i.e. Sv = 10 log10(sv)]. All calculations were performed using linear units.
Metrics for description of the vertical distribution of pelagic biomass, as estimated by acoustic backscatter.
| Quantity | Metric | Symbol | Formula | Unit | Reference |
|---|---|---|---|---|---|
| Density | Mean volume-backscattering strength | Sv |  | dB re 1 m−1 | 1 |
| Abundance | Area-backscattering strength | Sa |  | dB re 1 m2 m−2 | 1 |
| Location | Centre of mass | CM |  | M | 2, 3 |
| Dispersion | Inertia | I |  | m−2 | 2, 3 |
| Occupied area | Proportion occupied | Pocc |  | – | 2, 3 |
| Evenness | Equivalent area | EA |  | m | 2, 3 |
| Aggregation | Index of aggregation | IA |  | m−1 | 2, 3 |
| Layer structure | Number of layers | Nlayers | See explanation in text | Layers | – |
| Quantity | Metric | Symbol | Formula | Unit | Reference |
|---|---|---|---|---|---|
| Density | Mean volume-backscattering strength | Sv | | dB re 1 m−1 | 1 |
| Abundance | Area-backscattering strength | Sa | | dB re 1 m2 m−2 | 1 |
| Location | Centre of mass | CM | | M | 2, 3 |
| Dispersion | Inertia | I | | m−2 | 2, 3 |
| Occupied area | Proportion occupied | Pocc | | – | 2, 3 |
| Evenness | Equivalent area | EA | | m | 2, 3 |
| Aggregation | Index of aggregation | IA | | m−1 | 2, 3 |
| Layer structure | Number of layers | Nlayers | See explanation in text | Layers | – |
In all formulae, z represents the depth, and sv(z) the volume-backscattering coefficient at depth z. H is the total water column depth. The numbers in the reference column refer to: 1, MacLennan et al. (2002); 2, Bez and Rivoirard (2001); 3, Woillez et al. (2007).
Metrics for description of the vertical distribution of pelagic biomass, as estimated by acoustic backscatter.
| Quantity | Metric | Symbol | Formula | Unit | Reference |
|---|---|---|---|---|---|
| Density | Mean volume-backscattering strength | Sv | | dB re 1 m−1 | 1 |
| Abundance | Area-backscattering strength | Sa | | dB re 1 m2 m−2 | 1 |
| Location | Centre of mass | CM | | M | 2, 3 |
| Dispersion | Inertia | I | | m−2 | 2, 3 |
| Occupied area | Proportion occupied | Pocc | | – | 2, 3 |
| Evenness | Equivalent area | EA | | m | 2, 3 |
| Aggregation | Index of aggregation | IA | | m−1 | 2, 3 |
| Layer structure | Number of layers | Nlayers | See explanation in text | Layers | – |
| Quantity | Metric | Symbol | Formula | Unit | Reference |
|---|---|---|---|---|---|
| Density | Mean volume-backscattering strength | Sv | | dB re 1 m−1 | 1 |
| Abundance | Area-backscattering strength | Sa | | dB re 1 m2 m−2 | 1 |
| Location | Centre of mass | CM | | M | 2, 3 |
| Dispersion | Inertia | I | | m−2 | 2, 3 |
| Occupied area | Proportion occupied | Pocc | | – | 2, 3 |
| Evenness | Equivalent area | EA | | m | 2, 3 |
| Aggregation | Index of aggregation | IA | | m−1 | 2, 3 |
| Layer structure | Number of layers | Nlayers | See explanation in text | Layers | – |
In all formulae, z represents the depth, and sv(z) the volume-backscattering coefficient at depth z. H is the total water column depth. The numbers in the reference column refer to: 1, MacLennan et al. (2002); 2, Bez and Rivoirard (2001); 3, Woillez et al. (2007).
Total backscatter, a proxy of total biomass in the water column, was measured using the area-backscattering strength Sa, the integral of volumetric backscatter (sv) over the entire water column expressed as a decibel value (MacLennan et al., 2002). Mean density was measured by the mean volume-backscattering strength Sv. Occupancy was calculated as the proportion (Pocc) of the water column with Sv above −90 dB re 1 m2 m−2. Mean location was measured using the centre of “mass” CM, the average of all depths sampled weighted by their sv values. Bez and Rivoirard's (2001) inertia I measures dispersion or spread as the sum of squared distances from the centre of mass, weighted by the sv at each distance and normalized by the total sa. Evenness was measured using the equivalent area EA, calculated as the squared integral of sv over depth (i.e. 
) divided by the depth integral of 
. Values represent the area that would be occupied if all datacells contained the mean density (Woillez et al., 2007). This quantity can alternatively be expressed as its reciprocal, the index of aggregation IA, which is high when small areas are much denser than the rest of the distribution. Metrics were calculated from the processed data using Python scripts with the SciPy module for efficient array operations (Jones et al., 2001-2011). Applying these spatial statistics through the water column for each ping yielded a set of time-series that summarized each quantity of interest over the entire deployment. Figure 2 shows a representative echogram over 2 d, with the corresponding metric time-series. The metrics track diel vertical migration (DVM), including the formation and dissolution of layers, as well as the passage of aggregations near the surface.
Echogram for 2–4 March 2009, with the corresponding metrics plotted below. From top to bottom, metrics are mean volume-backscattering strength (Sv), area-backscattering strength (Sa), centre of mass (CM), inertia (I), proportion occupied (Pocc), equivalent area (EA), index of aggregation (IA) and the number of layers (Nlayers). Refer to methods and Table 1 for calculation and detail.
The final metric in the suite, the number of scattering layers, used an image-analysis approach. We defined layers as local scattering maxima in the vertical direction and used the slope of backscattering intensity with depth to identify them. The echogram was first smoothed by convolution with a Gaussian kernel to eliminate small-scale variability. The first and second derivatives of Sv were then estimated in the vertical direction on the smoothed echogram by taking the first and second differences along the vertical axis. Layers were defined as locations where the absolute value of the first derivative was close to zero, and the second derivative was negative. The first of these conditions selected areas where Sv did not change much with depth, meaning the tops of ridges (i.e. layers) and the bottoms of the troughs between them. The second condition, a negative second derivative, ensured that only areas at the tops of ridges would be selected. Before counting the layers, two additional filtering operations were performed. The first was a convolution with a median filter, replacing each cell's value with the median of all cells in its neighbourhood. This operation ensured that pixels meeting the two slope-based conditions would not be counted as layers if they were isolated in “non-layer” regions. The second was a convolution with a dilation filter, replacing each pixel with the maximum value in its neighbourhood. This operation eliminated isolated “non-layer” pixels within “layer” regions. It was then possible to count the number of transitions between “layer” and “non-layer” in each ping. This number, divided by 2, yielded the number of layers in the water column.
The number of layers detected depends in part on the parameter values used to define the layers: the acoustic threshold, the slope threshold defining where the derivative is close to zero, and the widths of the three filtering windows. To select the appropriate values for these parameters, we visually inspected “type” echograms from different times of the year. When two people compared the number of layers identified by eye, the greatest difference between their counts was one layer, and always fewer than ten layers were identified visually. Initial values were selected for the algorithm's parameters by starting with the smallest possible filter widths and incrementally widening them until the algorithm appeared to give results agreeing with the number of layers identified by the eyes. Based on this heuristic approach, initial parameter values of 19 × 19 pixels (9.5 m × 95 s) were selected for the standard deviation of the Gaussian filter, 0.2 dB m−1 for the slope threshold, and 19 × 19 pixels for the median and dilation filters.
To test the response of the layer detection algorithm to different combinations of input parameters, we conducted a Monte Carlo-based global sensitivity analysis. Input parameters were randomly selected using a Latin hypercube sampling (LHS) design (McKay et al., 1979). In an LHS design, the possible range of each parameter is divided into n strata, sized according to the probability density, and one value is selected randomly from each of these intervals. In each iteration, one of these previously selected values is drawn randomly for each parameter, making the number of parameter sets tested in the Monte Carlo run equal to that of strata n. LHS sampling has the advantage of ensuring that all portions of each input parameter's range are tested simultaneously, without requiring the testing of all parameter permutations. The complete coverage of the multidimensional parameter space is not guaranteed in a single run, so multiple runs need to be conducted. Rose (1982) found that an LHS design with 200 runs produced stable results; this number is used here. Each parameter's range was divided into ten strata of equal size, with ranges spanning the practical width of each parameter's values. The three convolution filters could have widths between 1 and 31 pixels, representing smoothing operations ranging from no smoothing to the approximate width of the layers identified visually. The slope threshold, controlling which areas of the echogram were designated as “flat” (and hence candidates for inclusion in layers), could take values from 0 to 0.15 dB m−1. Based on informal testing, values above 0.15 dB m−1 included regions of the echogram that were sloped. The range of the acoustic threshold was set from −100 to −65 dB.
To ensure the applicability of the results to the entire dataset, we divided the acoustic record into 30 equal temporal strata and randomly selected one data file from each stratum to repeat the Monte Carlo runs. Each file contained acoustic data from a period of ∼3 h. Within each file, 30 pings (i.e. individual profiles of the water column) were selected randomly in which to count layers after the application of the detection algorithm. Testing data from all periods of the deployment permitted comparison of the variability attributable to parameter adjustments with variability attributable to changes in the layer structure over time. All other metrics underwent Monte Carlo sensitivity testing analogous to that used on the layer detection algorithm, although there was only one parameter to adjust, the acoustic threshold, making an LHC design unnecessary.
Pattern description
For computational convenience and to ensure even spacing in time, all metric series were averaged into 1-h bins. A principal component analysis (PCA) assessed the degree of colinearity among metrics. As metrics do not share the same units, they were transformed to have zero mean and unit variance before conducting the PCA. All metrics displayed an annual trend, which was removed by fitting sinusoids with 1-year and 0.5-year periods to the data using ordinary least squares. All further analyses were performed on the residuals from these fits.
Fourier and wavelet transforms were used to characterize the variability in the detrended metric series as a function of temporal scale. Wavelet transforms decompose a series across time and frequency simultaneously, allowing the analysis of series that are non-stationary or contain frequency components at some times and not others. This is accomplished by convolving the series with a localized waveform (i.e. the wavelet), which is stretched or compressed to different scales, analogous to periods in Fourier analysis. The product is a two-dimensional array, mapping variability in the series as a function of time and scale. Wavelets offer the advantage over Fourier transforms of being able to locate transient frequency components in time. They also offer advantages over the short-time Fourier transform, in that their resolution is not fixed (Kaiser, 1994), allowing good time resolution for high-frequency events and good frequency resolution for low-frequency events.
As both Fourier and wavelet transforms require gap-free series with evenly spaced observations, a gap-filling scheme was used to fill periods of missing data. We used the stochastic interpolation method described in Percival et al. (2008). This algorithm uses the autocorrelation structure of the observed data to interpolate the expected value of missing observations between points on either side of the gap, then adds a stochastic noise component with the same autocorrelation structure and frequency spectrum as the rest of the series. The simulated data in the gap can be regarded as a realization of one section of the stochastic process that produced the rest of the time-series. The gap-filled series can then be used to calculate any subsequent statistics, such as the Fourier or wavelet spectrum, requiring continuous observations. By generating many realizations of the stochastic component, calculating the desired statistic each time, then taking the mean of all calculated statistics, confidence in the results is increased. This method uses the best linear predictor of the missing values, given the observations on either side of the gap and the autocorrelation structure of the series. By preserving the periodic and stochastic characteristics of the rest of the series, it reduces bias in estimates of the Fourier and wavelet spectra.
To carry out this gap-filling procedure, an autocovariance function (ACVF) for the series is necessary. Percival et al. (2008) used theoretical ACVFs calculated from autoregressive (AR) and fractionally differenced (FD) time-series models fitted to their time-series. Based on preliminary examination, the data contained both periodic behaviour and long-range autocorrelation, meaning that reasonably low-order AR and FD models would be poor fits to the data. Rather than try to build a complicated time-series model, we took a non-parametric approach and used the sample ACVF. In all, 500 realizations of the gap-filling process were generated, with the Fourier and wavelet spectrum of each realization calculated as described below. The mean of these 500 spectra was used as the final value for each transform.
Power spectra were calculated for the gap-filled series using a fast Fourier transform (Cooley and Tukey, 1965), normalized by a factor of N−1σ−2 (where N is the length of the series and σ2 is its variance). The final spectrum, averaged over all 500 gap-filling realizations, was smoothed with a modified Daniell kernel of width 15 (Daniell, 1946). Regression curves were fit to the unsmoothed final spectrum to model the dependence of variability on frequency.
Wavelet spectra of the metric series were calculated using continuous wavelet transforms (CWTs), following procedures described in Torrence and Compo (1998). The CWT is a highly redundant (i.e. non-orthogonal) transform, although the redundancy gives it fine resolution in time and scale (Torrence and Compo, 1998; Cornish et al., 2006). There are a number of continuous wavelet functions available, each with different properties. We used the Morlet wavelet with a frequency of 6, which has been used in ecological analyses (e.g. Ménard et al., 2006; Keitt, 2008) and offers a good compromise between time and frequency resolution (Torrence and Compo, 1998). The wavelet spectra were tested for significance at the 0.95 level against the theoretical spectrum of a first-order AR process fit to each series (Torrence and Compo, 1998).
Results
Metric sensitivity and characteristics
The layer detection algorithm proved relatively insensitive to the choice of parameters near the initial values selected visually (Figure 3). The average number of layers detected decreased gradually from 10.6 to 3.1 as the acoustic threshold was raised from −100 to −65 dB. This result reflected the exclusion of weaker scattering layers below a higher threshold. As the slope threshold was raised from 0 to 0.08 dB m−1, the mean number of layers detected increased from 0.9 to 8.0, then remained relatively flat through the maximum threshold tested, 0.3 dB m−1. The algorithm was most sensitive to the widths of the convolution filters used, in particular the Gaussian smoothing filter and the median filter. The mean number of layers increased from 0.4 to 10.5 as the standard deviation of the Gaussian filter was increased from 1 to 8 pixels, then decreased to 5.6 as the filter was widened to the maximum value tested, a standard deviation of 30 pixels. When the median filter was at its narrowest, the mean number of layers detected was 16.9, decreasing as the filtering window was widened. Once the median filter was wider than 13 pixels, the mean number of layers remained relatively constant between 4.7 and 6.8. The number of layers was less sensitive to the dilation filter, decreasing approximately linearly from 8 to 3 as the filter was widened from 1 to 33 pixels. For all three filters, narrower widths gave less predictable behaviour. At their lowest values, some layer counts as high as 50 or 60 were observed, far outside the range of values identified visually. At wider widths (>10–12 pixels), the coefficients of variation (CV) for all three were relatively steady, with values <1. Looking at the number of layers across all parameter combinations as a function of the time of year, there was a clear seasonal pattern in mean number, matching the pattern identified visually (not shown).
Sensitivity of the layer detection algorithm to the choice of parameters. Box-and-whisker plots show mean, interquartile, and 95th percentiles of number of layers counted at each parameter value (left axis). Black points show the CV at each parameter value (right axis). Subplots show the response of the algorithm with different values of (a) acoustic threshold, (b) width of Gaussian smoothing filter, (c) zero-slope threshold, (d) width of median smoothing filter, (e) width of dilation filter, and (f) data file (i.e. data subset) on which the algorithm was tested.
The remaining metrics displayed different degrees of sensitivity to the choice of acoustic threshold (Figure 4). The values of CM and I were insensitive to the threshold, although their CVs increased and decreased, respectively, as the threshold was raised. Values of EA and Pocc both decreased as the threshold was raised. Total water column Sa decreased very slightly when the threshold was raised above −70 dB. At the same time, mean Sv increased, reflecting the exclusion of less-dense areas with low backscatter values. This masking of low-density areas was also reflected in the IA which increased when the threshold was raised above −70 dB.
Sensitivity of metrics to the choice of acoustic threshold. Box-and-whisker plots show mean, interquartile, and 95th percentiles of metric values (left axis). Black points show the CV at each acoustic threshold (right axis). (a) Mean volume-backscattering strength (Sv), (b) area-backscattering strength (Sa), (c) centre of mass (CM), (d) inertia (I), (e) proportion occupied (Pocc), (f) equivalent area (EA), (g) index of aggregation (IA).
The PCA of the metric series showed a moderate degree of collinearity (Figure 5). The first component accounted for 62.0% of total variance, with the second and third components accounting for 17.2 and 8.6%. When summed, the first four principal components accounted for 95% of the variance. Pocc, Nlayers, and I were all correlated and aligned with the first component, and CM was aligned negatively with the first component. The IA and EA were anticorrelated with each other, with loadings orthogonal to those of Sa and Sv, which were strongly collinear.
Principal component analysis of metric time-series. (a) Biplot, showing the values of the first two principal components (points) and the loadings of the metrics along these two components (labelled). (b) Same as (a), but showing the second and the third principal components. The numbers in parenthesis on the axis labels show the proportion of variance explained by the corresponding component.
Pattern description
Collectively, the eight metrics captured a seasonal cycle in the vertical density distribution of Monterey Bay's pelagic fauna (Figure 6). Total integrated backscatter reached a minimum (Sa < −55 dB) in May of both years of the deployment, with the water column largely empty (Pocc < 0.1). Biomass was concentrated near the surface (CM> −200 m). By late summer, Sa had increased to approximately −40 dB because of the appearance of a large, non-migrating scattering layer between −400 and −600 m. The presence of this deep layer lowered the CM to −500 m and raised the occupancy to 0.7. Over autumn and winter, the deep layer gradually thinned, causing the occupied proportion of the water column to decline and the CM to rise. Total backscatter and mean density (Sa and Sv) remained high during that period, owing to the presence of a dense scattering layer near the surface. By spring 2010, the vertical distribution was similar to that seen in the previous year, with Sa 14 dB below its winter peak and concentrated high in the water column (i.e. CM near −100 m). From its lowest point in May 2010, Sa climbed rapidly through June and July, coinciding with the reappearance of the deep scattering layer. There was typically a distinct boundary at the bottom of the lowest scattering layer, below which only low densities of individual targets were observed. This boundary moved upwards in late winter and spring to a depth of ∼400 m before moving down again in summer, reaching depths of >700 m in early September 2009 and late August 2010.
Metrics time-series for entire deployment (27 February 2009 to 18 August 2010). The grey lines show values binned at 1-min intervals, and the black line shows daily averages. (a) Mean volume-backscattering strength (Sv), (b) area-backscattering strength (Sa), (c) centre of mass (CM), (d) inertia (I), (e) proportion occupied (Pocc), (f) equivalent area (EA), (g) index of aggregation (IA), and (h) number of layers (Nlayers). Refer to methods and Table 1 for calculation and detail.
This seasonal cycle accounted for large proportions of the variability in all metrics. Regressions of metric values on annual and twice-annual sinusoids (Table 2) were all significant (p << 0.001). Adjusted R2 values indicated that these sinusoids accounted for between 24.6% (IA) and 87.7% (Pocc) of total variability. All other series had R2 values between 0.474 and 0.688. Even after de-seasonalizing the series, they remained positively autocorrelated with lags between 4 and 27 d.
Summary of seasonal models fitted to metric time-series, with seasonal cycles estimated as the sum of sinusoids, with periods of 1 and 1.5 years.
| Metric | Source | MS | F4, 9390 | p-value | R2 | Maximum | Minimum |
|---|---|---|---|---|---|---|---|
| Sv | Model | 7.787 × 103 | 2.693 × 104 | ≈0 | 0.534 | −70.8 | −65.9 |
| Error | 2.89 | 10 September | 5 May | ||||
| Sa | Model | 1.467 × 104 | 4.597 × 103 | ≈0 | 0.662 | −37.5 | −44.2 |
| Error | 3.19 | 9 September | 5 May | ||||
| CM | Model | 1.681 × 107 | 5.167 × 103 | ≈0 | 0.687 | −154.8 | −421.5 |
| Error | 3.253 × 103 | 7 May | 14 September | ||||
| I | Model | 8.299 × 010 | 2.118 × 103 | ≈0 | 0.471 | 33 581.6 | 17 222.6 |
| Error | 3.918 × 107 | 11 December | 11 April | ||||
| Pocc | Model | 24.25 | 1.683 × 104 | ≈0 | 0.878 | 0.796 | 0.519 |
| Error | 1.440 × 10−3 | 9 September | 4 May | ||||
| EA | Model | 1.627 × 107 | 2.976 × 103 | ≈0 | 0.559 | 408.8 | 145.5 |
| Error | 5.469 × 103 | 3 September | 1 May | ||||
| IA | Model | 1.651 × 10−2 | 7.690 × 102 | ≈0 | 0.246 | 0.0099 | 0.0024 |
| Error | 2.147 × 10−5 | 19 April | 17 August | ||||
| Nlayers | Model | 7.489 × 103 | 3.632 × 103 | ≈0 | 0.607 | 8.57 | 3.56 |
| Error | 2.062 | 17 September | 21 April |
| Metric | Source | MS | F4, 9390 | p-value | R2 | Maximum | Minimum |
|---|---|---|---|---|---|---|---|
| Sv | Model | 7.787 × 103 | 2.693 × 104 | ≈0 | 0.534 | −70.8 | −65.9 |
| Error | 2.89 | 10 September | 5 May | ||||
| Sa | Model | 1.467 × 104 | 4.597 × 103 | ≈0 | 0.662 | −37.5 | −44.2 |
| Error | 3.19 | 9 September | 5 May | ||||
| CM | Model | 1.681 × 107 | 5.167 × 103 | ≈0 | 0.687 | −154.8 | −421.5 |
| Error | 3.253 × 103 | 7 May | 14 September | ||||
| I | Model | 8.299 × 010 | 2.118 × 103 | ≈0 | 0.471 | 33 581.6 | 17 222.6 |
| Error | 3.918 × 107 | 11 December | 11 April | ||||
| Pocc | Model | 24.25 | 1.683 × 104 | ≈0 | 0.878 | 0.796 | 0.519 |
| Error | 1.440 × 10−3 | 9 September | 4 May | ||||
| EA | Model | 1.627 × 107 | 2.976 × 103 | ≈0 | 0.559 | 408.8 | 145.5 |
| Error | 5.469 × 103 | 3 September | 1 May | ||||
| IA | Model | 1.651 × 10−2 | 7.690 × 102 | ≈0 | 0.246 | 0.0099 | 0.0024 |
| Error | 2.147 × 10−5 | 19 April | 17 August | ||||
| Nlayers | Model | 7.489 × 103 | 3.632 × 103 | ≈0 | 0.607 | 8.57 | 3.56 |
| Error | 2.062 | 17 September | 21 April |
The table shows the analysis of variance and estimated maximum and minimum values of the fitted seasonal model, with the corresponding dates.
Summary of seasonal models fitted to metric time-series, with seasonal cycles estimated as the sum of sinusoids, with periods of 1 and 1.5 years.
| Metric | Source | MS | F4, 9390 | p-value | R2 | Maximum | Minimum |
|---|---|---|---|---|---|---|---|
| Sv | Model | 7.787 × 103 | 2.693 × 104 | ≈0 | 0.534 | −70.8 | −65.9 |
| Error | 2.89 | 10 September | 5 May | ||||
| Sa | Model | 1.467 × 104 | 4.597 × 103 | ≈0 | 0.662 | −37.5 | −44.2 |
| Error | 3.19 | 9 September | 5 May | ||||
| CM | Model | 1.681 × 107 | 5.167 × 103 | ≈0 | 0.687 | −154.8 | −421.5 |
| Error | 3.253 × 103 | 7 May | 14 September | ||||
| I | Model | 8.299 × 010 | 2.118 × 103 | ≈0 | 0.471 | 33 581.6 | 17 222.6 |
| Error | 3.918 × 107 | 11 December | 11 April | ||||
| Pocc | Model | 24.25 | 1.683 × 104 | ≈0 | 0.878 | 0.796 | 0.519 |
| Error | 1.440 × 10−3 | 9 September | 4 May | ||||
| EA | Model | 1.627 × 107 | 2.976 × 103 | ≈0 | 0.559 | 408.8 | 145.5 |
| Error | 5.469 × 103 | 3 September | 1 May | ||||
| IA | Model | 1.651 × 10−2 | 7.690 × 102 | ≈0 | 0.246 | 0.0099 | 0.0024 |
| Error | 2.147 × 10−5 | 19 April | 17 August | ||||
| Nlayers | Model | 7.489 × 103 | 3.632 × 103 | ≈0 | 0.607 | 8.57 | 3.56 |
| Error | 2.062 | 17 September | 21 April |
| Metric | Source | MS | F4, 9390 | p-value | R2 | Maximum | Minimum |
|---|---|---|---|---|---|---|---|
| Sv | Model | 7.787 × 103 | 2.693 × 104 | ≈0 | 0.534 | −70.8 | −65.9 |
| Error | 2.89 | 10 September | 5 May | ||||
| Sa | Model | 1.467 × 104 | 4.597 × 103 | ≈0 | 0.662 | −37.5 | −44.2 |
| Error | 3.19 | 9 September | 5 May | ||||
| CM | Model | 1.681 × 107 | 5.167 × 103 | ≈0 | 0.687 | −154.8 | −421.5 |
| Error | 3.253 × 103 | 7 May | 14 September | ||||
| I | Model | 8.299 × 010 | 2.118 × 103 | ≈0 | 0.471 | 33 581.6 | 17 222.6 |
| Error | 3.918 × 107 | 11 December | 11 April | ||||
| Pocc | Model | 24.25 | 1.683 × 104 | ≈0 | 0.878 | 0.796 | 0.519 |
| Error | 1.440 × 10−3 | 9 September | 4 May | ||||
| EA | Model | 1.627 × 107 | 2.976 × 103 | ≈0 | 0.559 | 408.8 | 145.5 |
| Error | 5.469 × 103 | 3 September | 1 May | ||||
| IA | Model | 1.651 × 10−2 | 7.690 × 102 | ≈0 | 0.246 | 0.0099 | 0.0024 |
| Error | 2.147 × 10−5 | 19 April | 17 August | ||||
| Nlayers | Model | 7.489 × 103 | 3.632 × 103 | ≈0 | 0.607 | 8.57 | 3.56 |
| Error | 2.062 | 17 September | 21 April |
The table shows the analysis of variance and estimated maximum and minimum values of the fitted seasonal model, with the corresponding dates.
Variability was also present on subannual time-scales. The Fourier and wavelet power spectra of the deseasonalized series had most of their energy concentrated at low frequencies (i.e. long periods/scales), with the energy decreasing roughly as a power of the frequency. Power-law fits on the Fourier spectra had exponent values ranging from −1.732 to −0.8399, with R2 values between 0.18 and 0.31 (Figure 7). This energy was not distributed equally through time, as shown by the wavelet spectra (Figure 8). Significant energy at scales between 60 and 130 d was present during spring 2010 in all metrics except the index of aggregation, which reflected large changes in animal density distributions from March through June. In April 2010, there was a strong frequency component at the 30–40 d scale in Sv, Sa, and Pocc. Local peaks in the wavelet power were also present for most metrics at scales between 4 and 30 d, but they rarely attained statistical significance.
Fourier power spectra of all metrics, with the variance in each time-series plotted as a function of frequency (log axes). Subplots are labelled as in Figure 6, the grey points show the raw spectra, and the black lines show the spectrum after smoothing with a 15-point modified Daniell kernel. Regression lines show power-law distributions fitted to the raw spectra, with their exponents (i.e. the slope of the line on the log–log plot) and R2 values displayed. Vertical lines are plotted at 24- and 12-h periods (left and right lines, respectively).
Wavelet transforms of the metric time-series, plotting the variability in the series as a function of the time and temporal scale. Shading indicates the magnitude of the wavelet power, normalized by 1/σ2. The thick black contour line encloses regions significant at the 95% level, and dotted lines indicate the “cones of influence”: wavelet coefficients outside (i.e. below) these lines may be reduced in magnitude by edge effects.
DVM of scattering layers was observed nearly year-round, although it was not uniform through time. The power spectra of all metrics except IA had significant peaks at the 24-h period, as well as near the 12-h period in the CM, Pocc, and I spectra (Figure 7). Significant peaks near 12 and 24 h were also seen in the wavelet spectra of these series, particularly during February and March in the CM and I. The strong daily oscillation in the metrics at those times is attributed to the presence of several layers migrating to the surface from depths >300 m. The non-migrating component of the deep layer was mostly absent then. Significant variability was also observed in the wavelet spectra at time-scales <12 h, mainly during winter for CM, EA, I, Pocc, and Nlayers, and in Sa and Sv during their spring minima. During February and March of both years, dense aggregations were present near the surface during daylight, typically passing through the beam in several minutes and registering as sharp spikes in Sa, Sv, CM, and IA.
Discussion
The seasonal cycle in density distributions of animals in the water column coincided with the seasonal oceanographic cycle of Monterey Bay and the California Current. The spring backscattering minima were just before the typical peak of upwelling in June, when cool, nutrient-rich water shoals close to the surface, with subsequent peaks in chlorophyll and primary production (Pennington and Chavez, 2000). Within 3 months of its minimum, backscatter had increased an order of magnitude, with a much greater proportion of the water column occupied. This agrees with previous observations in Monterey Bay that peak zooplankton biomass lags primary production by 2–3 months (Robison et al., 1998; Croll et al., 2005). This lag is due either to the growth of these populations following the increase in food availability or to aggregation at the upwelling front, which is offshore early in the season but collapses shorewards in late summer (Abbott and Barksdale, 1991; Robison et al., 1998; Croll et al., 2005). Interestingly, the total biomass (represented by total acoustic backscatter Sa) remained high through winter and peaked in January, when primary production is typically close to zero (Pennington and Chavez, 2000).
Seasonal movements of the boundary between the mostly occupied upper water column and the mostly empty lower water column also appeared consistent with the known seasonal cycle. This boundary was between −400 and −700 m, consistent with the upper edge of the oxygen minimum zone (OMZ) in Monterey Bay. The OMZ, usually defined by concentrations of dissolved oxygen <0.5 mg ml−1, is located between ∼−500 and −1000 m and moves upwards during the upwelling season (Lynn et al., 1982; Robison et al., 2010). Physiological constraints imposed by low oxygen concentrations contribute to the structure of ecological communities in the meso- and bathypelagic depth zones (Robison, 2004). Our observations showed that the depth range of sound-scattering animals is restricted by oxygen availability.
Characteristics of the frequency spectra suggested that fluid turbulence influences the temporal variability of pelagic fauna. Spectrogram slopes of CM, I, Pocc, Sv, and Sa were all close to the theoretical value of −5/3 for diffusive turbulent mixing (Platt and Denman, 1975), and power-law distributions provided good fits to these empirical spectra. The observations indicated that the total abundance, mean vertical location, and dispersion of animals in the water column varied as the velocities of passively drifting particles. In contrast, the variance of IA and Nlayers decayed more slowly with increasing frequency. The temporal variability of these metrics, both measuring aggregation, was not consistent with passive-particle variability, confirming that midwater animal assemblages are mobile and aggregated independent of fluid motions. This finding is consistent with other observations of the distributions of mobile aquatic organisms.
Considerable variability was also present at short time-scales. There were significant frequency components with a 24-h period in all metrics except IA corresponding to the diel cycle. Significant variability at a near-12-h period was also present in the CM and Pocc indices, likely because of movements of animals by tidal currents. Diel variability was more pronounced in location and dispersion metrics than in those indexing total abundance and density. This contrast is expected if DVM is primarily vertical. The fact that Sa and Sv displayed energy peaks at the 24-h period indicates some cyclic variability in total abundance and density at that temporal scale. Variability at that scale could be attributable to changes in acoustic target strength as animals change orientation during DVM (Foote and Ona, 1987; Simard and Sourisseau, 2009) or to a horizontal component in the diel migration cycle (e.g. Benoit-Bird and Au, 2004). It could also be due to the migration of animals into the top 10 m of the water column, which were excluded from analysis.
Because there was no direct sampling, we cannot positively identify the organisms observed from DEIMOS. Monterey Bay has been sampled for a number of decades, and we can infer species compositions from the published literature. Barham (1956), in the first acoustic study of the Bay's ecology, identified the main constituents of acoustic backscatter layers as myctophid fish (Diaphus theta and Lampanyctus leucopsarus), and to a lesser extent, mesopelagic shrimp (Sergestes similis) and Pacific krill (Euphausia pacifica). Later studies, including visual observations from ROVs, have shown that layers are of mixed composition (Robison, 2004; Robison et al., 2010) and that as much as one-quarter of the pelagic biomass is gelatinous (Robison, 2004). Particularly notable are organisms such as the siphonophore Nanomia bijuga, a common gelatinous predator that is also a strong acoustic target as a result of its gas-filled pneumatophore (Barham, 1963; Robison et al., 1998). ROV observations show that Nanomia undergoes a seasonal cycle in abundance lagging primary productivity by ∼3 months, with peak densities near 500 m deep (Robison et al., 1998). Single targets in the region below the backscattering layers are likely macrourids and other fish (Yeh and Drazen, 2011). Future studies using stationary acoustics would benefit from direct sampling and additional acoustic frequencies, allowing better discrimination of different types of scattering organism (Horne, 2000).
Vertically applied spatial metrics other than Sv and Sa are relatively uncommon in fisheries acoustics. A typical approach would subdivide the water column into layers of interest and track the mean density or total backscatter within each of these layers. This was the approach taken by Flagg et al. (1994) and Cochrane et al. (1994) to add detail to the basic backscattering statistics. Burgos and Horne (2008) rigorously investigated a wide range of metrics to characterize walleye pollock (Theragra chalcogramma) aggregations in the Bering Sea, arriving at a similar set to that used in this paper.
Although the metrics do not measure compositional variability (i.e. changes in the abundance of different taxa), they do provide a powerful, high-scope view of aggregate variability (i.e. changes in the overall abundance or density; cf. Micheli et al., 1999). As demonstrated by the PCA, considerable collinearity exists between metrics, as would be expected, given the mathematical relationships between them (Bez and Rivoirard, 2001; Woillez et al., 2007). As metrics were calculated over the entire water column, they integrated distinct features observed in the acoustic record. For example, the signature of diel variability attributable to vertical migration in the upper water column was damped in the metrics by the presence of the deep non-migratory layer. The CM and I indices both assume a single centre, but the vertical density distribution was typically multimodal, with several distinct layers. Depending on the process or the region of interest, the metrics could be constrained to specific regions in the water column.
The suite of metrics presented here effectively captures the vertical density distribution, rendering a large dataset into a tractable form, and provides a flexible set of tools to quantify variability in the ecosystem. They are capable of measuring periodic variability (e.g. DVM), long-term shifts (e.g. interranual changes), and transient events (e.g. dense aggregations). These phenomena are important when characterizing ecosystems and could be used to detect potential changes associated with climate trends. Metrics used in this study could be extended and applied in two or three spatial dimensions as succinct descriptors of density distributions through time. The metrics are well-suited to pattern description in high-scope dataseries: they are objective, simple to calculate, retain full resolution in the indexing dimension (i.e. time), and provide a description of the vertical density distribution with a minimum of assumptions. Other applications potentially include unsupervised monitoring of a datastream or as preliminary indicators of data quality.
Long-term, continuous sonar records can be used to address ecological questions. DVM is an obvious process to investigate with stationary acoustics, but appropriately located acoustic instruments could also be used to observe the timing of horizontal migrations, either on short (diel) time-scales (Benoit-Bird and Au, 2004) or on longer scales for more mobile organisms (Kaltenberg et al., 2010). When combined with oceanographic observations, stationary acoustics facilitates the investigation of biophysical coupling through the water column across temporal scales. Stationary acoustics could also be used to complement ship-based surveys for resource management as methods for temporally based stock assessment are developed (cf. Brierley et al., 2006). Stationary acoustics offer a powerful tool to monitor an ecosystem at multiple trophic levels (Koslow, 2009). The efficiency and efficacy of monitoring is greatly improved with a suite of objective metrics such as that presented here.
Acknowledgements
We thank Kongsberg Maritime for the loan of the echosounder used in this study, and Dick Kreisberg for his role in the design and assembly of DEIMOS. We also thank MBARI for hosting DEIMOS at the MARS observatory and, in particular, the crew of the RV “Point Lobos” and ROV “Ventana”. We greatly appreciated the helpful comments on the draft manuscript provided by editor Verena Trenkel and three anonymous reviewers. Funding was provided by the University of Washington School of Aquatic and Fishery Sciences.
References
Author notes
Handling editor: Verena Trenkel
